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COMPUTER PROGRAMS FOR PRESSURIZATION (RAMP) AND PRESSURIZED 
EXPULSION FROM A CRYOGENIC LIQUID PROPELLANT TANK 

by Philip A. Masters 
Lewis Research Center 

SUMMARY 

An analysis to predict the pressurant gas requirements for the discharge of cryo- 
genic liquid propellants from storage tanks is presented, along with an algorithm and two 
computer programs. One program deals with the pressurization (ramp) phase of bring- 
ing the propellant tank up to its operating pressure. The other program deals with the 
expulsion of the liquid at a uniform pressure. The method of analysis involves a numer- 
ical solution of the temperature and velocity functions for the tank ullage at a discrete 
set of points in time and space. The input requirements of the program are the initial 
ullage conditions, the initial temperature and pressure of the pressurant gas, and the 
time for the expulsion or the ramp. Computations are performed which determine the 
heat transfer between the ullage gas and the tank wall. Heat transfer to the liquid inter- 
face and to the hardware components may be included in the analysis. The program out- 
put includes predictions of the mass of pressurant required, the total energy transfer, 
and the wall and ullage temperatures. 

The analysis, the algorithm, a complete description of input and output, and the 
FORTRAN IV program listings are presented in this report. Sample cases are included 
to illustrate the use of the programs. 


INTRODUCTION 

Planning for space vehicle and mission requirements necessitates continuing opti- 
mization of propellant tank pressurization systems. This optimization is realized in an 
accurate determination of pressurant requirements for any given set of operating param- 
eters, such as tank pressure, inlet gas temperature, liquid outflow rate, and tank size. 
This knowlec^e will allow the design of a system that carries only the weight (pressurant 
gas and associated tankage) necessary to accomplish the mission. 



The analysis in reference 1, which is the basis for this report, provides a selected 
set of simplifying assumptions to predict the quantity of pressiirant gas required during 
the pressurized discharge of liquid hydrogen. Evaluations of pressurant gas injectors 
based on this work are made in reference 2. These evaluations deal with small-scale 
tests made in a 0. 82 -cubic -meter (29-ft ) cylindrical tank. The reference 1 analysis 
proved to be adequate (within ±10 percent) in predicting the pressvmant gas requirements 
even though two of its major assumptions were shown to be invalid - namely, no heat 
transfer to the liqmd surface and no mass transfer. 

The purpose of this report is to revise and extend the analysis of reference 1 to in- 
clude the energy transfer occurring at the gas -liquid interface in tanks of any arbitrary 
symmetric shape. The analysis was also modified and extended to cover the initial pres- 
surization (ramp) period. The major limiting assumptions still remaining from the orig- 
inal analysis are one -dimensional flow and the exclusion of mass transfer. 

The algorithm developed in the analysis employs numerical solutions to the gas tem- 
perature and velocity functions for the tank ullage at a discrete set of points (in space and 
time) called net points. Two separate computer programs, coded in FORTRAN IV, are 
given (each based on the use of a single -component pressurant): one for the pressurized 
expulsion of a cryogenic (single component) liquid from an axisymmetric storage tank, 
and the second for the pressurization (ramp) period which may precede an expulsion. 

The programs are used independently, although they may be coupled together to predict 
pressvirant gas requirements as well as the ullage gas temperature distribution and the 
adjacent tank wall temperature distribution for the entire pressurization cycle. 

Although prior knowledge of the operating conditions for a fluid system is needed for 
an analysis, the use of the analytical programs does not necessitate that ejqperimental 
data derived from prototype systems be available. 

The validity of the analysis presented herein has been verified in references 3 to 6 
for 1. 5- and 4 -meter (5- and 13 -ft) diameter spherical tanks using either gaseous hy- 
drogen or helium as the pressurant over a range of inlet temperatures, tank pressures, 
and outflow rates. The predicted pressurant requirements for the expulsion period, for 
all cases, were 7.0 to 12.4 percent of the measured (experimental) values. The pre- 
dictions for the ramp phase were ±0. 5 to ±16. 0 percent of the measured values. 

The analysis, the algorithm, a complete description of input and output, and the 
FORTRAN IV listing are given in this report. Sample cases are included to illustrate 
the use of the programs. 


ANALYSIS 

The general method of analysis is developed in reference 1 and briefly summarized 
here. The general analytical model is shown in figure 1. It is based on an application of 
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the first law of thermodynamics to the fluids and materials within a cryogenic -liquid 
storage tank which has the pressure history shown in figure 2. The first law is repre- 
sented by a general energy eqviation which is coupled with a transformed, one- 
dimensional equation of continuity. Substitutions involving the equation of state and the 
equation of heat transfer are made to effect the transformation. 

In the analytical model shown in figure 1, the ullage volume in a partially filled liq- 
uid propellant tank is divided into a set of volume elements. This is done by establishing 
a series of equally spaced net points along the vertical axis of the tank to the liquid inter- 
face. In each volume element, convective heat transfer from (1) gas to wall, (2) gas to 
liquid, and (3) gas to internal hardware is treated in a steady-state manner by an equa- 
tion which expresses gas temperature as a function of velocity and position. 

The analysis treats the ullage as a single -component gas. For the ramp program, 
the ullage pressure is made functionally dependent on a discrete set of pressure -against - 
time coordinates which is submitted as part of the input data. The initial set of net 
points remains fixed throughout the ramp period and the hold period which may follow. 
For an expulsion, liquid propellant is expelled, as pressurant gas is added, to maintain 
the pressure within set limits in the ullage. A new net point is added to the ullage for 
each time increment necessary to advance the liquid interface as the liquid is discharged. 

To determine a imique solution to the velocity and temperature functions at each net 
point, initial and boundary conditions are specified. These conditions are described fully 
in the section INPUT -OUTPUT REQUIREMENTS. 

The form of the energy equation used in the analysis in reference 1 and modified to 
account for both arbitrary symmetric tank shapes and internal tank heat sinks may be 
written as 


ar 

at 


^^c (T^ - T) 

rMPCp 



V aT ^ ^^^1 3P ^ ^ 

ax jMPCp at VpCp VpCp 


( 1 ) 


where 

T^ temperature of tank wall 
% heat -transfer rate to internal hardware 

Qj^ heat -transfer rate at liquid interface 

(All symbols are defined in appendix A. ) 

By a substitution involving the equation of state, the one-dimensional equation of 
continuity is transformed into a functional relation involving the velocity of the ullage gas 
as a function of temperature and position. For the one -dimensional expression of 
continuity. 
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± (pVA) + 1 (pA) = 0 
3x 9t 


( 2 ) 


2 

The substitution A = Trr is made, where r is the position radius at location x along 
the vertical axis. The expression for density from the equation of state p = MP/ZRT is 
also substituted: 



= 0 


(3) 


The following velocity equation is obtained after performing the partial differentiation 
and after rearrangii^ terms: 


3x 




1 ^ 

p at 


9r 

r 0x 


(4) 


Each of the bracketed terms may be simplified by differentiating the equation of state 
while holding the pressure constant and again while holding the temperature constant. 


Z 


1 


Z 


2 



When the ejq)ressions for 


Zj and Z 2 are substituted into equation (4), 


dx ZT\at 9x/ zp at r ax 


(5) 

( 6 ) 


(7) 


where the final term on the right is the contribution of the tank curvature. 
The heat -transfer equation at a point in the tank wall can be written 


aT 


w 




(T - T^) + 


iW 




( 8 ) 


where q^ is the rate of heat addition per imit area to the tank wall from outside the 
tank. 
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The equations contained in the analysis are too complex for a closed -form solution. 
The numerical method used here and described in the algorithm is brot^ht about by ap- 
proximating the differential equations by algebraic equations. For example, the pre- 
ceding equation is approximated by 


%,j ' . ^ 9 — \ (t! -t^ j)+ ( - 0 

At Vw^W^wZ ’ Vw^w^w/ 

J J 

where the prime refers to the new value of the variable. Solving the equation for T^ j 
gives 



h At 
c 




% 


h At 


1 + 


V^w^w^w/ 


( 10 ) 


where the quantities marked with the asterisk may be evaluated at the beginning or the 
end of the time interval. The subsequent algorithms will proceed to show the transfor- 
mation resvilting when T^ . from equation (10) is substituted into the algebraic approxi- 
mation of equation (1). The solution of the resulting transformed equation is coupled to 
the algebraic approximation of equation (7), and the solution techniques are described. 

There are two parts to the algorithm in each computer program. The preliminary 
part deals with input data and problem definition. The main computation part is con- 
cerned with the determination of the amount of pressurant gas added dxiring the expulsion 
or pressurization period. 

The region of the distance -time plane in which the solution is carried out is defined 
by a set of net points equally spaced along the vertical axis. The entire set of net points 
is based on an assumed value of a specified number in the initial ullage. From experi- 
ence, an assigned number of 5 to 6 net points per decimeter (15 to 20 net points/ft) of 
ullage is satisfactory for nearly all situations and generally provides the desired accu- 
racy for a reasonable computer execution time. Too many net points in the initial ullage 
may result in using all the available storage for the variables before the expulsion is 
completed. 
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CALCULATION PROCEDURE 


Expulsion Algorithm 


A step-by-step description of the basic calculation procedure is given here. Steps 1 
to 6 refer to the program listing shown in appendix B under PRELIMINARY COMPUTA- 
TION. For the solution to proceed, a set of boxmdary and initial conditions is required. 
These conditions are specified in the section INPUT -OUTPUT REQUIREMENTS. Steps 7 
to 17 describe the main computation and the results. Figure 3 gives a logic diagram for 
the expulsion algorithm. A listing of the expulsion program is given in appendix B. 

Step 1: make units conversion, make geometry calculations, and interpolate initial 
wall and gas temperatures and specific heats . - The program is structured so that, at 
the user’s option, the input data are printed out in SI units or in U. S. customary units. 
First -value parameters are set equal to zero, and values for the constants are computed. 
The space between the points is established by setting the initial number of net points for 
the designated initial ullage. 

Step lA begins with a logic statement which provides three basic options to the pro- 
gram user. For a cylindrical tank, the radius is specified. Under this option the pro- 
gram solution is based on a tank with hemispherical end sections joined to the cylindrical 
midsection, as would be specified in the input data. If the end sections are not hem- 
ispherical, the radius is set equal to zero. Input coordinate values are specified (see 
Group II data in section INPUT -OUTPUT REQUIREMENTS), and the program interpolates 
the radius for each volume element. The same procedure is followed for spheroidal 
tanks; however, for a spherical tank, the tank radius is set equal to -1 and the coordi- 
nates for each volume element are determined from the tank diameter specified in the 
input (group I data) . 

Figure 1 establishes the program model and the tank configuration on a set of co- 
ordinate axes. The total number of points selected for the ullage is part of the input 
data. Although each point is located at the center of its element, each variable asso- 
ciated with the point is representative of the entire element. The distance separating the 
points, once fixed at the start, remains constant for all points throughout the run time. 
The first boundary point is at the top of the tank, where x = 0. The other boundary point 
is at the liquid-vapor interface, and each boundary element is one-half the thickness of 
the other elements. For an expulsion the interface is advanced one point for each time 
step, which is determined by the propellant discharge rate. 


dt = dx 


A(N„) 



( 11 ) 
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For a spherical tank, the radius and flow area for each volume element at the 
point j are 



The tank weight for the configuration is approximated by the relation 


= 






( 12 ) 

(13) 


(14) 


The wall thickness is evaluated by the interpolation subroutine from the discrete 
values of wall thickness as a function of distance from the top of the tank. 

At the initial time (t = tQ) the gas temperatures and the wall temperatxires as well as 
the specific heat for each net point in the initial ullage are defined based on a linear inter- 
polation of the input data in step 1C. 

Step 2: compute initial values of heat -transfer coefficients. - Gas transport prop- 
erties at the mean of the gas and wall temperatures for each net point are computed and 
the free -convection correlation 


h„ = - n(GrPr)*” (15) 

L 

is employed to evaluate the heat -transfer coefficient at the point. Program options allow 
for input of multiplier n and e:^onent m. The effect of the flow length L is canceled 
when the default value of 1/3 is used for the exponent since the ler^th L is raised to the 
third power in the Grashof parameter. The pressurant gas properties are included for 
the subroutines where the computation is made; however, the coefficient h^ may be 
specified as a constant for an entire rvin. 

Step 3: compute initial values of q^, Q (internal hardware), P, and AP/At. - 
The outside-wall heating rate, the inside -hardware heating rate, the tank pressime, and 
its first time derivative are initialized. Typical values in references 3 to 6 are 

= 10 W/m^ (0. 0009 Btu/ft^-sec) and Q = -8. 7 W/m (-0.0025 Btu/sec -lineal ft) . A 
subsequent option under step 9 is provided in which the energy transferred to the internal 
tank hardware is computed. The hardware component temperatures are initialized here 
if the option is taken. If the option is not taken, positive input data values are negated in 
the program. 
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step 4: compute initial value of compressibility factor and then its derivatives which 
are needed in continuity equation . - The local compressibility factor for the gas temper- 
ature and tank pressure, as well as its derivatives, is defined. The equation of state for 
a real gas is commonly written in terms of a compressibility factor Z(P,T) 


1 _ ZRT 
p MP 

Upon differentiating and holding pressure constant. 


where the expression 



Differentiating again but with respect to P, holding T constant. 


where 



Z 



= Z 


2 


(16) 


(17) 


(18) 


(19) 


( 20 ) 


Step 5: compute initial values of local ullage gas velocities . - The parameters and 
temperatures evaluated in the previous steps are introduced into a substituted form of 
the continuity equation. This is obtained by substituting equation (1) into equation (7) and 
noting that Is zero at time tj = 0. The equation is first solved at the point adjacent 
to the interface, and the solution is continued from point to point to the top of the tank. 
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The differential equation (eq. (25), ref. 1) has been modified here to include spher- 
ical tanks and approximated at net points (Xj, tj) by 


■ 1 

Ax 2 


^2h„ZiR 

_c ^ 

rMPC 
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( 21 ) 


Step 6: find gas in ullage at time zero by a numerical integration over ullage density 
profile . - The mass of the volume elements in the initial ullage is totaled. Since the in- 
terface is located at the center of its volume element and the first net point is a boundary 
condition, by definition, only one-half the masses in these two elements are included in 
the total. 



( 22 ) 


where = f(T^> ■ 

Step 7: compute initial calculations . - All necessary input data required for the 
main part of the calculation are now available. The identification of key parameters, 
along with input and boundary conditions, is important when the output resvilts are 
reviewed. 

Step 8: find temperatures at new time. - When T^ j from equation (10) is substi- 
tuted for the value of T^ in equation (1), 



The equation is equation (Cl) in references 3 to 6, in which the quadratic is eval- 
uated for the gas temperature T1 beginning with the second net point (j = 2) at the 
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top and the value for t! ^ as the boundary value. The solution to this equation is 
repeated from point to point until all the ullage temperatures are computed. The value 
of Tj j is always the temperature evaluated at the previous point for the new time. The 
variable Tj is the temperature evaluated for the point at the previous time. 

The program computes the real positive root of equation (23) by use of the quadratic 
formula. The wall temperature T^ j is also computed under step 8 by equation (10), 
and the value required for Tj is obtained from the gas temperatvire quadratic (eq. (23)). 

Step 9: evaluate energy transfer to internal hardware . - A program option is pro- 
vided in which the energy transferred to internal tank hardware may be computed. The 
computed value for step 9 wovild default any input value. The new hardware temperatures 
are computed by using the relation 


T^ = h 




A„ — 3iLi . At + T 


(MCv)^ 




The heat transferred is then 



(24) 


(25) 


If four hardware components make up the hardware, the average heat -transfer rate is 



Step 10: find compressibility factors at new time . - In step 10 the compressibility 
factor and its derivatives are reevaluated based on the new gas temperatures computed 
in step 8. 

Step 11: find velocities at new time . - The new gas velocities are reevaluated by 
using the new temperature profile and compressibility factors. The finite difference 
form of equation (4) is used. Althovigh the ullage temperatures are computed starting at 
the top of the tank, the velocity equation is used to calculate the ullage gas velocity 
starting with the point Ny (in fig. 1) near the interface. The velocity at the interface 
N^, the boundary value, is determined from the propellant discharge rate which is part 
of the input data. 
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( 27 ) 



Step 12: find heat flow rate to wall and total heat added to wall; find gas flow rate 
and to tal gas added up to new time . - The heat transfer to the tank wall is computed by 
using the relation 


Q - 


Q = 





(28) 


For the elapsed time tj -tg, the wall temperature changes from j j 


k . - 27 t Ax jPw 


(29) 


In step 12, the mass of pressurant in the tank is dependent on the outflow time as 
well as on the pressurant gas temperature and pressure. The mass of gas in the tank at 
any time t dviring outflow is based on a summation of the mass in all the volume ele- 
ments in the ullage. The initial mass calculated at the start is subtracted, and the differ- 
ence represents the pressurant mass added up to time t. 



(30) 


An alternate method for determining the flow of pressvu-ant into the tank is provided. 
The GASCHK parameter measures the velocity of the pressiu’ant passir^ by an arbitrary 
area element near the top of the tank. With the gas density evaluated at the element, the 
mass flow into the tank for a given expulsion is summed over the entire expulsion time. 
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t=i,f ^ 

Step 13: find specific heats at new time . - The wall and gas specific heats from the 
top of the tank to the interface are evaluated at the new temperatures. 

Step 14: find heat -transfer coefficient at new time . - In step 2 the use of the free- 
convection correlation is described. The gas-to-wall heat-transfer coefficient is ob- 
tained from the same correlation in step 14. This correlation is used for the conditions 
outlined in references 1 to 7. Reference 8 verifies that the gas-to-wall heat -transfer co- 
efficient is in the free -convection regime for a helium expulsion experiment but reports 
that it is definitely within the forced-convection regime for oxygen test data. In step 14 
this heat -transfer coefficient has been extended to include two component effects, if de- 
sired; a forced -convection component, as well as a natural -convection component. The 
following equations represent the manner in which reference 8 relates both convection 
components: 

Free -convection correlation: 


= 0.00117 r^ 

(32) 

h„ = li n(GrPr)“ 

(33) 

^ L 


Forced-convection correlation: 



(34) 

k \A^M/ \ k / 


^g,W " *^c ^0® 

(35) 


The last term in the combined equation represents the diminishii^ effect of forced con- 
vection as the distance x increases from the pressurant distributor. 

The development of the gas -to -liquid heat -transfer coefficient is given in appendix B 
of references 3 to 6. It was shown that the conductance across the gas -liquid interface 
is similar in form to the empirical relation for free -convection flow across a horizontal 
surface. 

h L 

Nu = = 0. 14(GrPr)”^ (36) 

k 
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This eqmtion, with a value of m = 1/3, is used in step 14 to coincide with the references 
cited. The heat transferred to the liquid is given by 




Tl,s> 


(37) 


in which Tg is a representative temperature at the edge of the thermal boundary layer. 
For hydrc^en or helium pressurant over liquid hydrogen, Tg was determined ejq)er- 
imentally to average 1. 3 times the adiabatic compression temperature given by 


T 


ad 




,(y-l)/y 


(38) 


Although the magnitude of is computed from a AT representii^ a gradient across 
a gas -liquid interface, the energy transferred to the liquid is represented in the energy 
equation (1) as being uniformly derived from the entire ullage. 

Step 15: check to see if end of time has been reached . - The time increment is com- 
puted from the displacement of the interface to the succeeding net point. The time incre- 
ments are summed to give the elapsed time, which is compared to the time specified for 
the discharge. 

Step 16: end of time exceeded - interpolate conditions at end of time . - When the 
sum of the time steps has exceeded the specified discharge time , a back interpolation of 
the computed values is performed to comply with the end of time specified for the 
discharge. 

Step 17: write out results . - The subroutine WRITE 2 is called and the computed 
values are printed. 


Pressurization (Ramp) Algorithm 

A separate computer program determines the mass of pressurant added, as well as 
tank wall energy requirements, during the ramp and hold periods. The same equations 
that describe the expulsion period are also applicable for the ramp and hold periods. 

Even though experimental results (refs. 3 to 6) indicate relatively large amounts of mass 
transfer during this period, mass transfer and heat transfer to the liquid are not included 
in the analysis because of the added complexity of the transport process occurring at the 
interface. A logic diagram for the ramp algorithm is given in figure 4. 

The analysis computes the gas temperatures in the ullage at any time durii^ the 
pressure rise (and hold period) from the gas energy equation. The corresponding gas 
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velocities are computed from the equation of continuity. An iterative technique is used 
here which describes how convergence is achieved in the solution for velocity. 

The predicted mass of pressurant added is based on a summation of the mass of the 
volume elements by assuming a one-component ullage consisting of 100 percent pressur- 
ant at the end of the ramp and hold periods after subtractii^ the initial mass at the start. 
This is explicit for an autogenous pressurization system, for example, where pressur- 
izing gas is derived from an engine jacket bleedoff. For a nonautogenous system, 
where the pressurizing gas differs from the propellant, a two -component ullage mixture 
is encountered once ramping is initiated. For a two -component mixture, the contribution 
to the final pressvtre by the component gas in the initial ullage is evident from the heat of 
compression. On a molar basis, the programmed procedime for obtaining the mass ad- 
dition by the difference of the assumed sii^le -component pressurant in the initial ullage 
from that in the final ullage would still appear satisfactory. This assumption is verified 
in references 5 to 7. 

Steps 1 to 7 : These steps follow the procedure described in the section Expulsion 
Algorithim, with the exception that hardware component temperatures are not initial- 
ized under step 3 since this computational option is not provided. A program listii^ for 
the ramp program is provided in appendix C. 

Step 8: find temperature at new time, find miscellaneous quantities (Qjj and q^) at 
new time, initialize an estimate for velocity. - The inlet gas temperature representing 
the upper boundary temperature is defined for the time t > 0. Likewise, the value of the 
saturation temperature is defined for the time t > 0 and is made equal to the lower gas 
boundary temperature at the interface. The botindary wall temperature at the interface 
is made equal to the bulk propellant temperature. 

The option to compute heat transfer to the internal hardware is not included 
during the ramp, since this effect during the ramp period is generally small. However, 
a substituted value for this parameter as well as the outside wall parameter (heat leak 
rate) q^ may be made a part of the input data. 

The initial gas velocity distribution used in solving the temperature function at each 
net point for each new time t is obtained from the previous time as follows: 


AP 


1-2 




At 


1-2 


AP 


0-1 


At 


0-1 


(39) 


A solution to the ullage temperature function (expressed as a quadratic in the expul- 
sion analysis) under the conditions encountered diiring ramp proved to be extremely 


14 


difficult. When the initial wall temperature distribution in the ullage was greater than 
the gas temperature distribution, heat transfer from wall to gas resulted in negative ve- 
locities, which made it impossible to evaluate the function. A satisfactory method is 

achieved, however, when the velocity at the net point V. is eliminated from the fxmc- 

— ^ ] 

tion. Substituting the velocity fxmction for Vj into the quadratic results in the following 
cubic equation (appendix C of refs. 3 to 6): 


bjTf .1 


- 1 ) 


* \z /j * ix\z ). At 








(40) 


The main purpose of step 8 is the evaluation of the preceding eqiiation. Unlike the 
method employed in the expulsion program, the cubic equation is solved numerically by 
the Newton -Raphson method. 



where F(t!) represents the cubic fimction evaluated for the temperature T at the new 
time (’) for the point j. The term F*(T.') is its first derivative. The subscript i for 
T! represents the initial or previovisly evaluated temperature in the iterative process 
and^Tj^^ is the new value, which has converged closer to the true value (see fig. 5). 

Step 9: find heat -transfer coefficient at new time . - Step 9 computes the gas -to -wall 
heat-transfer (natural convection) coefficient. Steps 9 to 11 together allow for the deter- 
mination of the thermodynamic and transport properties as a function of P and T in the 
subroutines. 

Step 12: find velocities at new time . - The velocity equation used to calculate ullage 
gas velocity starting with the point (fig. 1) is the same equation employed in step 11 
of the expulsion algorithm. The boundary value of the velocity at the interface, point 
N^, is zero with no expulsion. Figure 5 shows the iteration scheme by which conver- 
gence is achieved in the solution of the gas energy and continuity equation. 
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The ullage gas velocity is calculated from point to point until the top of the tank is 
reached. The new velocities are used in the cubic equation in step 8 along with the pre- 
vious values of and the temperature distribution is redetermined. This process 

is repeated until the computed velocities are essentially unchanged from those computed 
in the previous iteration. When convergence is achieved over the entire ullage, the time 
is then advanced to tg and the process is continued as shown in figure 4. 

Step 13: find gas flow rate and total gas added up to new time . - Generally, the 
ramp rate is not accompanied by tank outflow. In step 13, the mass of pressurant in the 
tank is dependent on the pressurant gas temperature and pressure as a function of ramp 
time. The mass of gas in the tank at any time t during the ramp is based on a sum- 
mation of the mass in all the volume elements in the ullage. The initial mass calculated 
at the start is subtracted, and the difference represents the pressurant mass added up to 
time t 



t=f t=0 


Nf 

Ni 

S Pn "^^n 

“ ^ Pn ^^n 

n=l 

n=l 

t=f 

t=i 


(30) 


An alternate method for deter minii^ the flow of pressurant into the tank is provided. 
The GASCHK parameter establishes a velocity for the pressurant passii^ through the 
second area element from the top of the tank. With the known gas density at the cross 
section, the mass flow into the tank diming the ramp is summed over the entire ramp 
time. This method generally did not give good agreement with the computed mass de- 
rived from the integration of all the volume elements 


M = 



At 


(42) 


Step 14: find heat flow rate to wall and total heat added to wall . - This step follows 
the description for step 12 in the expulsion program. 

Step 15: check to see if end of time has been reached . - The time increments are 
summed to give the elapsed time, which is compared to the specified time for the ramp 
or hold pressure. In general, ramp has never been experienced with significant tank 
outflow. However, this possibility does exist, and the time step would then become a 
function of the discharge rate. Under this condition, the procedure outlined by this 
step becomes severely constrained. The time involved for the interface to advance to 
the next point is in many instances too great and results in a serious discontinuity. 
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Nevertheless, the option to include tank discharge durii^ ramp is provided and may be 
considered as an extension of the expulsion program with an iteration scheme which may 
provide some utility for a pressvirization system under study. 

Step 16: end of time exceeded - interpolate conditions at end of time; and step 17: 
write out results . - In essence, steps 16 and 17 follow the corresponding steps of the ex- 
pulsion program . 


SUBROUTINES 

The subroutines, which are common to both the pressxirization and expulsion pro- 
grams, are listed in appendix D, 


Subroutine SINTS 

SINTS is a subroutine which allows the user the option to program in SI units. This 
subroutine converts the input data in SI units to the units required by the program. 


Subroutine SPHEAT 

The ’’specific heat” subroutine SPHEAT provides 15 storage locations representii^ 
a set of 15 discrete values of specific heat as a function of temperature (°R) for the ves- 
sel wall, as well as a set of 15 values of specific heat (Btu/lbm-'^R) as a function of tem- 
perature (°R) for the pressurant gas. These data may be substituted to fit the problem 
definition. The program interpolates values between the points. 


Subroutine WRITE 1 

Subroutine WRITE 1 contains the statements which print out the input data. 

Subroutine WRITE 2 

Subroutine WRITE 2 contains the statements which print out the problem solution 
(the calculated pressurant gas and temperature distributions) at any time t during the 
expulsion and/or the ramp, following the call statement. 
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Subroutine COMPRS 


COMPRS is a subroutine which calculates the compressibility factor from a 20x17 
matrix for hydrogen (Z(T,P)). For another gas pressurant the table may be replaced, 
or a value of Z(T,P) = 1. 0 or Z(T,P) = Constant may be substituted. 


Subroutine HCOEFF 

Subroutine HCOEFF calculates heat -transfer coefficients for either of two 
pressurant gases (hydrogen or helium) at a specific pressxme. Transport property data 
for each gas are keyed to the molecvilar weight and may be replaced for other gases. 


Subroutine INTERP 

Subroutine INTERP performs a straight-line interpolation between two discrete 
points. 


INPUT -OUTPUT REQUIREMENTS 
Input 

The program input consists of a three -card description of the problem, althov^h any 
or all of the cards may be blank. Any specific information may be entered in columns 2 
to 80, but column 1 is left blank. The succeeding groups of data are defined in the expul- 
sion program listing. 

Group I data . - Group I data include the parameters (PARAMS) that relate to tank 
configuration and input options. These data are entered in NAMELIST FORM, 

Group n data . - Group II data are also entered in NAMELIST FORM, which uses 
names in place of FORMAT numbers in the read (INPNTS) and write (OTPNTS) state- 
ments. The INPNTS series is a sequence of coded statements (words) giving the number 
of pairs of values read in as data which define the initial and boundary conditions for the 
ramp or expulsion programs. The flexibility of the NAMELIST FORM is demonstrated 
in the write statement (OTPNTS), which verifies the number of pairs of values that are 
specified. The number of pairs in each set printed under OTPNTS must agree with the 
number of paired values in each set submitted under the NAMELIST TABLES that follow. 
If a sii^le pair is specified under TABLES at time t = 0 (i. e. , TGAS = 400. , TIMEl 
= 0. ,), the value of the variable (TGAS) is maintained constant for any subsequent time. 
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The expulsion program provides a maximum of 15 data entries under TABLES. The 
ramp program has the capability of utilizing only the first 11 data entries defined in the 
following list. The data for MRAD are necessary if the tank radius (as a function of dis- 
tance from the top) is to be interpolated. 


Code 


Number of pairs of data 


MTGAS 

MPDATA 

MFLOW 

MTSAT 

MTBULK 

MQOUT 

MQIN 

MTWIN 

MTGIN 

MRAD 

MTICK 

MTBAK 

MTCU 

MTSS 

MTAL 


inlet gas temperature (K; °R) as fvmction of time (sec) 

2 2 

tank pressure (MN/m ; Ibf/ft ) as function of time (sec) 
volume flow rate (m /sec; ft /sec) as function of time (sec) 
saturation temperature (K; °R) as function of time (sec) 
liquid -propellant bulk temperatvire (K; °R) as function of time (sec) 
outside-wall heating rate (J/m -sec; Btu/ft -sec) as function of time (sec) 
inside -hardware heating rate (J/m-sec; Btu/lineal ft-sec) as function of 
time (sec) 

initial wall temperature (K; °R) as function of axial distance from top of 
tank (m; ft) 

initial gas temperature (K; °R) as function of axial distance from top of 
tank (m; ft) 

tank radius (m; ft) as function of axial distance from top of tank (m; ft) 
tank thickness (m; ft) as function of axial distance from top of tank (m; ft) 
initial temperature of phenolic internal hardware (K; °R) as function of 
axial distance from top of tank (m; ft) 
initial temperature of copper internal hardware (K; °R) as function of 
axial distance from top of tank (m; ft) 
initial temperature of stainless steel (TP304) internal hardware (K; °R) as 
function of axial distance from top of tank (m; ft) 
initial temperature of aluminum internal hardware (K; °R) as function of 
axial distance from top of tank (m; ft) 


Group in data. - The term group m data refers to data that are entered directly into 
the appropriate subroutine. 


Output 

The output from a successfully executed case is written by the subroutine WRITE 2 
after a printout of the initial problem conditions by subroutine WRITE 1. This output, 
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shown for a sample case in appendixes B and C , consists of a block of computed wall 
temperatures and gas temperatures as a function of distance above the liquid interface. 

The block data are preceded by two lines of data. The top line gives the time 
(TIME) during the ramp or expulsion for which the temperatures are evaluated. Asso- 
ciated with the time, the GAS SUPPLIED, HEAT TO WALL, and INLET VEL (pressurant 
velocity) are included on the first line. For an expulsion, the second line of output in- 
cludes the heat transfer to the liquid (propellant). Both programs provide the CHECK 
ON GAS computation, and the GAS FLOW or rate of pressvirant into the tank. 

Under group I data, the ENDTIM value specifies the ullage condition at the point in 
time for the termination. For the ramp program the ENDRMP value represents the 
ullage condition at the end of the press\ire rise period, and ENDTIM represents the con- 
dition at the end of the hold period (fig. 2). 

The OUTPUT parameter specifies the number of time steps for a ’’demand” require- 
ment of the computational results. 

Two techniques are used to determine the mass of pressurant required for any given 
time t. The primary method determines GAS SUPPLIED by summit^ the product of gas 
density and volume for all the volume elements for each net time. From this summed 
value, the initial ullage mass is subtracted. The alternate method, CHECK ON GAS, 
establishes a velocity for a given cross section near the tank inlet. With the gas density 
known at the cross section, the mass flow rate into the tank is summed over the entire 
time. 


TYPICAL ERROR MESSAGES AND THEIR CAUSES 
Ramp Program 

Initial ullage velocities are vmstable . - Since the wall and gas temperature distribu- 
tions are specified in the initial ullage, the values for gas specific heat, heat -transfer 
coefficient, and compressibility factor are defined. The initial velocity at each net point 
is computed for time t = 0 and P = PDATA(l), where PDATA(l) is the initial pressure 
throughout the ullage for a series of discrete time -dependent values. This initial ve- 
locity distribution is printed out by the program. A negative value for a velocity indi- 
cates an incipient instability, and the program terminates itself because of a major over- 
flow. The initial problem condition must be altered by selecting a greater number of net 
points (smaller net space between points) or by adjusting the initial slope of the pressure 
rise curve (fig. 1). 

Solution to cubic equation requires more than 10 iterations . - Sometimes during the 
solution for the true ullage temperature, the Newton -Raphson programmed convergence 
is not achieved in 10 iterations. These occurrences are counted and when the number is 
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greater than 25, the program reduces the time step and starts the problem again. If the 
solution is not achieved by this procedure, the pressure rise for the interval is too great, 
or the instability may be attributed to a large difference between the upper -boundary gas 
temperature and the inlet pressurant temperature. 

One or more negative gas velocities are computed after time zero . - If, dxiring the 
ramp, one or more negative velocities are computed, the program reduces the time step 
and starts the problem again. It is programmed to do this three times and then prints 
out all the gas temperatures in the ullage. 

A computed temperature is greater than the boimdary temperature . - Sometimes 
during the course of the solution of the cubic equation, there will be three real and un- 
equal roots. Unfortunately, a distinction cannot be made as to the correct root. When the 
roots are near the boundary temperature or the inlet temperature, occasionally the con- 
verged value may be greater than the value of the inlet temperature. In this predicament, 
the root can be made equal to the value of the boundary temperature without interrupting 
the procedxire. 

Error greater than acceptable (program cannot converge on true ullage tem- 
perature) . - Path B in the ramp block diagram, shown in figure 4, should not be con- 
strued as an alternate program path. It indicates a sequence of printed data values and 
their relative occiarrence during the iterations. 

If after 32 iterations, any of the ullage velocities differ from the value computed in 
the previous iteration by more than 1/2 percent, this condition is indicated by a printout 
of the velocities. 

After 40 iterations, the program may still not converge on the true ullage temper - 
atxires. In this circumstance, the deviation of the velocities from those computed in the 
previous iteration may be greater than 10 percent. This situation is indicated by a print- 
out of the ramp time, ullage temperature, and wall temperature distribution. 

These data blocks are printed during steep ramp rates or when extremely nonlinear 
portions of the ramp curve (pressure rise as a function of time) are encountered. The 
ramp time parameters PMl, PM2, PM3, and PM4 should be examined and the time step 
selected accordingly. However, when the time steps are made exceedii^ly small, ex- 
cessive computer run times are encovintered. Time steps of 0. 2 second or less, selected 
for a large percentage of the total ramp time, may consume more than 1/4 hour of com- 
puter time for a 40- to 50 -second ramp even though the ullage volume may be only 5 per- 
cent of the tank volume. 


Expulsion Program 

For the expulsion program, error messages have not been found necessary. In gen- 
eral, if the program is terminated, the reason often appears obvious upon examining the 
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input data format. If the input data format is correct, a program completion is assured 
when the initial gas temperature boundaries are equivalent to the liquid interface con- 
dition and the pressurant temperature at the start. The wall temperature at the liquid 
interface should approximate the propellant bulk temperature. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 15, 1974, 

502-24. 
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APPENDIX A 


SYMBOLS 


Engineering FORTRAN Units in program 
symbol name 

A AREA 

AAY 

BEE 

TEST - 


Description 

Tank cross section normal to vertical axis 

3 2 

For cubic equation y +Py +qy+r=0, 
AAY = (3q - P^)/3. 

BEE = (2P^ - 9pq + 27r)/27 

Evaluate BEE^/4 + AAY^/27 to determine 
the nature of the roots. 


b G8 



BQUAD 



CP 

Btu/abm-°R) 

^v’^W 

CPW 

Btu/(Ibm-°R) 

C 

G1 

CQUAD 

°R 



D 

DIA 

ft 


DISC 



d G6 ft/sec 

Gr 

ERRP 


Coefficient of first term of gas temperature 
function, eq. (40), represented by cubic 
equation in ramp program 

Second term of quadratic temperature func- 
tion in expulsion program 

Specific heat at constant pressure 


Specific heat at constant volume, C^; 
specific heat of wall material, 


a 


aw 



At q^ 


Final expression of quadratic temperature 
function in expulsion program 


Spherical or cylindrical diameter; vertical 
diameter when tank is a spheroid 


Discriminant of quadratic temperature 
function 


^2 Ax/ p' - P \ Ax 
Z At\ P’ / Z At 

Grashof number, L^p^g/3 AT/p^ 

A percent difference in velocity, computed 
in step 12 (ramp) , from value of previous 
iteration 
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Engineering 

symbol 

FORTRAN 

name 

Units in program 

Description 

g 


ft/sec^ 

Gravity acceleration 


H 

Btu/(ft^-sec-°R) 

Convective heat -transfer 
coefficient 

— 

HOPE 


Temperature function, cubic in 
propellant temperature TP 

— 

ISTAR 


Integer representing count of 
iteration 


KSTAR 


Counting integer - specifies the 
number of times the numerical 
solution for a gas temperature 
could not converge on the true 
value in 10 iterations 


KLAMP 


Integer counter equal to or less 
than 3 . The ullage temper - 
ature computations are re- 
peated so that temperature - 
dependent parameters are 
evaluated close to converged 
gas temperatures. 


LOOM 


Counting integer - if computed 
value for a gas velocity is less 
than zero, time step is 
reduced. 

J 

XJAY 

(ft-lbf)/Btu 

mechanical equivalent of heat 

L 

XL 

ft 

Flow length 

1 

TICK 

ft 

Wall thickness 

M 

GSRATE 

Ibm/sec 

Pressimant flow rate addition 

M 

XMOLEC 

lbm/(lbm-mol) 

Molecular weight of pressvirant 
gas 

AM 

GAS (GASB - GSTART) 

Ibm 

Amount of gas (pressurant) added 
by subtracting initial ullage 
gas from GASB 
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Engineering FORTRAN Units in program Description 

symbol name 

— GASCHK Ibm Mass of pressurant added by 

alternate calculation 

My GASB Ibm Summed value of ullage gas over all vol- 

ume elements for time into ramp or 
expulsion 

m HEXP Grashof, Prandtl exponent 

N N,NP Ntmiber of volume segments in tank, 

NP, refers to next time step 

N. to Nj J Summii^ index, i = initial, f = final 

NjtoN 2 Particular volume segments 

Nu Nusselt number, h L/k 

^>^0 PjPHOLD Ibf/ft^ Tank pressure; initial pressure 

AP PP - P Ibf/ft^ Differential pressure 

Pr - Prandtl number, C ju A 

PHE Coefficient of second term (square in 

TP) of dibic equation in ramp program 

Q QjCQTR Btu Heat transfer to wall ; heat transfer to 

liquid interface 

Q,Qy DQjQTR Btu/sec Heat -transfer rate to tank wall; heat- 

transfer rate to liquid 

AQl, AQ2 DQl, DQ2 Btu Heat transfer to hardware 

components 1 and 2 

q^ QOUT Btu/ (ft^ -sec) Heat -transfer rate to wall from outside 

tank 

R R (ft-lbf)/(°R)(lbm-mol) Gas constant 

ROOT 1 °R Real positive root of temperatmre 

quadratic in expulsion program 

Re Reynolds number, LVP/p 

r RAD ft Tank radius 

RR Final term of cubic equation in ramp 

program 
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Engineering 

symbol 

FORTRAN 

name 

Units in program 

Description 

T 

TP,TWP 

°R 

Ullage temperature; tank wall temperature 

T 

L,S 

TP(N+1) 

°R 

Temperature of the satvirated propellant 

— 

TPP 

°R 

Assigned temperature equal to inlet 
pressurant temperature 

AT 

TDIFF 

°R 

Differential temperature 

T 

^6 

TADD 

°R 

Temperature at edge of thermal boundary 
layer 

t 

TIME 

sec 

Time into ramp, hold period, or expulsion 

At 

DT 

sec 

Time increment 

— 

UU 


Coefficient of third term (TP term) of gas 
temperature function (ramp program) 

V 


ft^ 

Ullage volume 

V 

V 

ft/sec 

Gas velocity associated with a specific net 
point 

AV 


ftS 

Volume increment 

— 

VP 


Gas velocity associated with a specific net 
point at previous iteration 


VHOLD 

ft/sec 

Gas velocity passing volume element near 
top of tank, used in alternate method of 
calculatii^ pressurant going into tank 

^n 

XN 


Number of net points in ullage 

X 

X 

ft 

Coordinate in direction of tank axis 

Ax 

DX 

ft 

Distance between net points 

Z 

Z 


Compressibility factor 

^1 

Z1 



^2 

Z2 


“-’(81 
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Description 


Engineering FORTRAN Units in program 
symbol name 


a 


0 

% 

Y 

M 

P 

(Ji) 


D3 °R 


h At 

i+_e 

Pw^W^W 



BETA 1/°R 


Vise lbm/(ft-sec) 

RHO Ibm/ft^ 

D4 


Coefficient of thermal expansion 

Dimensional decay coefficient of ullage 
forced heat -transfer correlation 

Specific -heat ratio 

Viscosity 

Density 

^ ^ RZQl \ At 

V ^ Mirr^Ax Mirr^Ax / ^p^ 
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APPENDIX B 


LISTING OF EXPULSION PROGRAM 


C THIS IS THE EXPULSION PROGRAM 

C 

C 

c 

DIMENSION 

1 TGAS(30)*TIMEl(30)tP0ATA(30) vTlHEZISO) t 

2 ELOH(30)«TIHE3I30I .TSAT(30> tTlHE4(30)» 

3 TBULK(30)*TIME5(30)*OOUTD(30)tTlME6(30l* 

4 CINO(30)*TIME7(30>«ZAV(150)* 

5 THlNO<3O>.OISTll3O>.TGINOI3O)tOIST2(30>. 

6 X(150)*T( 150>.TPU50I.TmI130)*TNPU50) •VIISOI. 

7 CPI 150>«CPU(150)*H(150itZil50) «Z1I150) »Z2( 150) * 

8 RADI 150) «AREA( 150) •YRADI 150) *XRAO(150) >OROX( 150) 

C 

DIMENSION 

1 C5( 150)*C9( 150)«VTICK(50)*XXT(50) .TICK (150) 

2 ,TWBK(30)«01ST7(30)»TkiCU(30).0IST8(30).TSTN(30)*01ST9(30)» 

3 TWAL(30).0IS11(30).CPBK(150).TWB( 150) »TB( 150) •CPCU(150)« 

A THC( 150).TC(150).CPST(150).TST(150) .TZU50) .CPAL(i50) * 

5 TALS4 150) .TBCI 150) «TTT ( 150) .PRAM( 150) 

C 

COMMON 

1 X.T.TP.TW.TWP.V.CP.CPM.TGAS.TIHEl.MTGAS.POATA.TlMEZ.MPOATAt 

2 fLOW*TIME3«MFLOH»TSAT.TIME4.MTSAT*TBULK.TIME5»MTBULK.OOUTD. 

3 TIME6.MOOUT.OINO.T1ME7.MOIN.N.NP.XN.ULLAGE. RADIUS. 00. 

4 SPWGT.XMOLEC.GSTART.HCONST. TIME. GAS. O.GASCHK.GSRATE.OT. 

5 /CONST. H.HMULT.HEXP.YRAO.XRAO.MRAO.TNKMT.COTR.QIN.UNUMS. 

6 mtmin.mtgin.mtick.mtbak.mtcu.mtss.htal. 

7 TMIN0.01ST1.TG1N0.01ST2.YTICK *XXT . TMBK. 01 ST7 . 

8 TUAL.0IS11.THCU.01ST8.TSTN.D1ST9 
C 

COMMON 

1 CARAO.AOIFU.SPNSS.OIA.RCONST.TAOO.AB.AIC. A3S.A5B.WTB.WT1C. 

2 NT3S.NT5B 
C 

C 

c 

C READ 3 CAROS OF PROBLEM OESCRIPTION AND NRITE OUT. THERE 

C MUST BE THREE CAROS USED. ALTHOUGH ANY OR ALL OF THEM MAY 

C BE BLANK. LEAVE THE FIRST COLUMN OF EACH CARO BLANK AND 

C ENTER ANY INFORMATION IN COLUMNS 2 TO 80. 

C 

1 WRITE (6.100) 

00 2 J«1.3 
READ (5.101) 

2 WRITE (6.101) 
f 

C 
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c 

c 

f 

c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


ncmenclature for input data 

XN« NUMBER OF NET POINTS AT TIME ZERO (MUST BE 3 OR MORE) 

NET POINTS AVAILABLE 150 (SEE DIMENSION STATEMENTS) 
REOUIREO INPUT GROUP I DATA 
OUTPUTf NUMBER OF TIME STEPS TAKEN BEFORE EACH OUTPUT 
ULLAGE* INITIAL ULLAGE HEIGHT* (CANNOT BE ZERO) 

RACIUS* TANK RADIUS 

IF RADIUS »0* THE RADIUS IS INTERPOLATED BASED ON THE TANK RADIUS 
(VS DISTANCE FROM TOP) DATA READ INTO THE PRCGRAM. 

IF RACIUS - -1* THE TANK ISA SPHERE (DIA * DIAMETER) 

SPNGT* TANK WALL SPECIFIC WEIGHT 
SPWSS* SPECIFIC WEIGHT OF TANK LID MATERIAL 
ENOTIM* TIME AT WHICH OUTFLOW ENDS* SECONDS 
XHGLEC* MOLECULAR WEIGHT OF PRESSURIZING GAS 

ZCONST* COMPRESSIBILITY FACTOR (1. FOR IDEAL GAS*BLANK FOR REAL) 
HCCNST* HEAT TRANSFER COEFF. (BLANK IF H IS TO BE CALCULATED) 

IF HCCNST IS BLANK H WILL BE COMPUTED FROM THE EQUATION 
H=HMULT ♦ COND/XL ♦ (GRASHOF ♦ PRANOTL) ♦♦HEXP 
HMULT* CONSTANT IN ABOVE EQUATION (0.13 IS USED IF LEFT BLANK) 

HEXP* CONSTANT IN ABOVE EQUATION (0.333 IS USED IF LEFT BLANK) 

DIA* DIAMETER WHEN THE TANK IS A SPHERE OR CYLINDER. 

SPECIFY DIAMETER ON VERTICAL AXIS WHEN TANK IS SPHEROID 
RCGNST* INITIAL HEIGHT USED IN CALCULATION OF H 

TR1.TR2 GOVERNS THE MODE OF TRANSFER BETWEEN PRESSURANT GAS AND 
TANK WALL. IF TR2 IS BLANK AND TR1=1. • HEAT TRANSFER 
IS BY FREE CONVECTION* IF TR2=^1.* TRl IS BLANK* HEAT 
TRANSFER IS BY FORCED CONVECTION 
TADC* A TEMPERATURE AT THE EDGE OF THE THERMAL BOUNDARY 

LAYER TO DETERMINE THE DRIVING POTENTIAL* (TADD-TSAT). 
FOR HYDROGEN PRESSURANT OVER LIQ H2* TAOO WAS DETERMINED 
EXPERIMENTALLY TO BE 1. 2-1.5 TIMES THE ADIABATIC 
COMPRESSION TEMPERATURE 
SEE NASA TN 0-5336* 53B7 
UNUMS SET GREATER THAN 0. FOR Si UNITS 

QQQQ^ IF QQQQ IS BLANK* PROGRAM ASSUMES NO INTERNAL HARDWARE 

UNLESS THE QIND VS TIME/ VALUES ARE GIVEN UNDER INPUT 
DATA. IF QQQQ*1* THEN SOME OR ALL PARAMETERS 
AB THRU WT5B ARE SPECIFIED 

A8* EFFECTIVE AREA PHENOLIC HARDWARE EXPOSED TO THE 

PRESSURANT GAS IN THE VOLUME ELEMENT 
AlC EFFECTIVE AREA OF COPPER HARDWARE EXPOSED TO THE 

PRESSURANT GAS IN THE VOLUME ELEMENT 
A3S* EFFECTIVE AREA OF THE 304 SS HARDWARE EXPOSED TO THE 

PRESSURANT GAS IN THE VOLUME ELEMENT 
A5B* EFFECTIVE AREA OF ALUMINUM HARDWARE EXPOSED TO THE 

PRESSURANT GAS IN THE VOLUME ELEMENT 
WTB* WEIGHT OF THE PHENOLIC HARDWARE IN THE VOLUME ELEMENT 

WTIC* WEIGHT OF COPPER HARDWARE IN THE VOLUME ELEMENT 

WT3S* WEIGHT OF 304SS HARDWARE IN THE VOLUME ELEMENT 

WT5B* WEIGHT OF ALUM. HARDWARE IN THE VOLUME ELEMENT 
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c OPTION 2 DATA 

QIC. PARAMETER USED IN CALCULATING GAS TO MaLL. 

FORCED CONVECTION HEAT TRANSFER 
AOIFU PRESSURANT DISTRIBUTOR AREA 

CARAO« CHARACTERISTIC RADIUS OF THE TANK WHEN FORCED CON- 
VECTION IS A MODE OF HEAT TRANSFER 
END CF OPTION II DATA 

BEXPO. MBEXPO PARAMETERS INVOLVING DECAY COEFFICIENTS 

SEE NASA TM X - 53165 


GROUP I DATA 

NAMELIST / PARAMS / XN «OUTPUT • ULLAGE *RADIUS » SPWGT* ENOT IMt XMOLECt 
lHCONSTfDlA«RCONST*ZCCNSTf HMULT«HEXP*AB*AlC«A3StA5B»WTB*MTlCt 
2hT3S.WT5B«SPMSS«0D0Uf TA0D*ClC*TRl«TR2tCARA0*A0IFUtCYLN.UNUMS 
REAO(5«PARAMS> 


C 

C 

C 

C 

c 


f. 

c 

c 

c 

C 


C 

c 

c 

c 


GROUP II CATA EXPULSION DATA- CONTINUED 


NAMELIST / INPNTS / MTGAS«MPDATA.HFLOW*MTSAT *MTBULKtNOOUT» MOINt 
lHThIN*MTGIN«MRAD«MTICK*MTBAK*MTCU*MTSS*MTAL 
READ! 5. INPNTS) 

NAMELIST / OTPNTS / MTGAS«MPDATA*MFLON *MTSAT «MTBULK»MOOUTtMOINt 
lMTbIN*MTGIN.MRAD*MTICK*MTBAK»MTCU*MTSS.MTAL 


WR ITE(6*0TPMS) 
NAMELIST / TABLES / 


1 

TGAS 


TIMEl 

9 

PDATA 

f 

TIME2 

f 

PLUM 

9 

TIME3 

9 

? 

TSAT 

♦ 

TIME4 

9 

TBULK 

9 

TIME5 

f 

OOUTO 

9 

TIME6 

9 

3 

UiNO 

• 

TIME? 

9 

TwlND 

9 

OISTl 

• 

TGINO 

9 

DIST2 

9 

4 

YRAO 


XRAO 

9 

YTICK 

9 

XXT 

f 

T^BK 

9 

OIST7 

9 

5 

TtoCU 

9 

OIST8 

9 

TSTN 

9 

0IST9 

9 

TWAL 

9 

DISH 



READ (5*TABLES) 


GROUP III PROPERTY DATA- SEE SUBROUTINE- 
NOTE - THE CATA BLOCKS IN THE SUBROUTINES MAY BE SUBSTITUTED FOR 
CIFFERENT GASES OR NALL MATERIAL CONSISTENT WITH 
HYDROGEN /HELIUM BEHAVIOR* 

1 - SUBROUTINE SPHEAT - MATERIAL AND GAS SPECIFIC HEATS 

NOTE — The SPECIFIC HEAT SUBROUTINE FOR HYDROGEN 
HAS NO PROVISION FOR PRESSURE VARIATION 
NOTE ~ THE WALL SPECIFIC HEAT DATA IN THE SUBROUTINE 
DOES NOT PROVIDE STORAGE FOR MORE THAN ONE MATERIAL 

2 - SUBROUTINE COMPRS - GAS COMPRESSIBILITY 

3 - SUBROUTINE HCOEFF - TRANSPORT PROPERTIES 


**«*«•«*«***#«* ««*** 
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c PMELIMINAFY COMPUTATION 

ITER1«0 

0»0. 

6ASChK«0* 

H2h=0.0 

OTR=0.0 

COTR»0.0 

TNKmTsO.O 

GSTART^O. 

OOlsO. 
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u on no 


c 

c 

301 

C 


4 

5 

C 

6 

199 


002^0* 

003^0. 

004=0. 

CALL WRITEl 

IftUNUMS .EO. 0. ) GO TO 301 

CALL SINTS 

CONTINUE 

XPI=3. 14159 

R»1545.4 

XJAV=778.2 

HCCNST=HC0NST/3600. 

N-XN 

NP»N+l 

0X=ULLA6E/(XN~1.0) 

IWR ITE=OUTPUT 


OIFU = S0RT(4.«ACIFU/XP1> 
C8=R/XM0LEC 
C3 = OX 4 C8 
Cl = C3/2. 

C2 = Cl/XJAY 
C4 » DX/C8 
C6 = 2,*C8 
C7 = C8/XJAY 


DO 3 J=2fl50 
XJT iCK=J-l 
X( J)*CX*XJTICK 

CALL INTERP ( YT ICK« XXT fHT ICK t X ( J ) * TIC K( J ) 1 
C5( J 1=1. 0/4 TICK U)*SPWGT) 

C9 4 J)=2.04XPI4>0X«T1CK(J)*SPMGT 
X( 1 ) = 0.0 

TICK! n=YTICK( 11 


STEP-IA- 

HAKE GEOMETRY CALCULATIONS 
IF (RACIUS) 8«4«6 
no 5 J=2*150 

CALL INTERP 4 YR AO* XRAO.HRAQ ,X 4 J> t R A0( J 1 1 

AREA4J>= XPI«RAC4 J)*«2 

KAC(11=RA0421 

AREA4 1 )=AREA421 

GO TO 12 

CONTINUE 
no 7 J»2.150 

IF4X4J1 - OIA/2. 1199*299. 299 
RACUl = S0RT40IA « X4J1 - X4J1««21 
AREA4J1 = XPI * RA04J14>«2 
GO TO 7 
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799 IFCXIJ) •GT. (CYLN * RAOILS)) GO TO 399 
RAC(J) » RADIUS 
AREAUl » XPl * RAO(J)*«2 
GO TO 7 

399 PLEG « X(J) - CVLN - RADIUS 

IPiPLEG .GE. RADIUS) GO TO 499 
RACUi - SORT! RAO IUS*«2 > PLEG**2) 

AREAIJ) s XPI * RAO(J)*4>2 

7 CCKTINUE 

499 RAC4 1 ) « RAC(2) 

AREA(l) » AREAI2) 

GO TO 12 

8 DO 10 J»2*1S0 

IE (X(J)>OIA) 9*9*11 

9 RAOIJ) - SORTIDIA « XU) > XU)**2) 

10 AREAU)aXPl«X(J)*(OIA-X( J) ) 

11 AREA! 1 )«AREA(2> 

RA0(l)sRA0(2) 

C 

12 CONTINUE 

DO 16 J«2*149 

13 ORCXlJ)«lRAOU>l)-RAOU))/tX( J-»l)-X( J) ) 

C 

C TANK WEIGHT DOES NOT INCLUDE WEIGHT OF LIU CP CONNECTOR AT TOP 

C 

15 TNKWT > TNKWT * C9I J)4>RA0(J)«S0RT( 1.0 0R0X(J)**2) 

16 CONTINUE 

17 0R0Xm«IRA0l2)-RA0I 1) )/X(2) 

0RCXI150) « 0R0XI149) 

0T*0XWAREA(N)/FL0WI1) 

CALL INTERP I FLOW«T I HE3* MFLOW * OT *FLO) 

OT^OXWAREAIN)/! FLOW! I ) -t-FLO )42« 0 
VP«FL0W(1)/AREAIN) 

TIME»CT 

C STEP-IB- 

C CCHPUTE INITIAL WALL AND GAS TEMPERATURES 

C 

00 19 J>1*N 

CALL INTERP ( TW INO.O I STl * HTW 1 N *X( J ) « TW U ) ) 

18 CALL INTERP I TGINO*OIST2* MTGIN *X( J ) . T( J) ) 

C 

C 

C STEP-IC- 

C COMPUTE INITIAL VALUES OF SI ECIFIC HEAT 

C 

19 CALL SPHEAT I T ( J ) *TW I J ) »CP U) * CPWl J) * XMOLEC) 

C 

C 

C STEP-2- 

C COMPUTE INITIAL VALUES OF HEAT TRANSFER COEFFICIENT 

C 

IF (HCONST) 1*22*20 
C 

20 00 21 J»l*150 

21 H(J)sHCCNST 
GO TO 31 
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c 

27 IF (HHULT) 1.23*24 

23 HMULTaO.13 
C 

24 IF (HEXP) 1.25.26 

25 HEXP*0.333 

f. 

26 00 30 J~l,N 

IF (XI J)-RC0NST) 27.27.26 

27 XLsOIA 
GO TO 29 

26 XLsOIA 

29 TTT(J)aThlJ) 

30 CALL HCOEFF ( T (J > «TTT ( J) . POATA ( I) . XL .HI J ) »PR AH ( J ) . ZCONST .HMULT .HEX 
IP.XHOLEC) 

31 CONTINUE 
C 

C STEP-3- 

C COMPUTE INITIAL VALUES UF UOUT. OIN. P AND OPOT 

C 

OOUTaCOUTCI II 
OIN a-CINOIll 
PapOATAI 1 ) 

PFCLOaP 

IF IMPCATA-1) 1*32.33 

32 PPaP 

60 TO 34 

33 CALL INTEPP I POATA. T IHE2.MP0AT A. T I HE . PP) 

34 OPDTal PP-PI/OT 

IFICOCO - 0* >144.144.42 

42 00 43 J*1*N 

CALL INTEPP ITWBK.OISTT.HTBAK.XI Jl .TWBI Jll 
CALL INTEPP (TMCU.OISTB.MTCU.XUl.TMCIJ) I 
CALL INTEPP (TSTN*CIST9*HTSS.XU).TST(JI ) 

43 CALL INTEPP I TMAL *01 SI l.HTAL « X I J I . TALS I J 1 1 

144 CONTINUE 

C 

C STEP-4- 

C COMPUTE INITIAL VALUES OF COMPRESSIBILITY FACTOR AND DERIVATIVES 

C 

IF iZCCNST) 1.37*35 

35 00 36 Jal.150 
7IJ|al.0 

21 1 J)al.O 

36 Z2(J)al.O 
ZhOLOal.O 
GO TO 39 

C 

3/ 00 3B Jal*N 

3B CALL CCMPRS ( T (J I .POATA 1 1 ) *2 1 J > *Z 1 1 J> * 22 (Jl . XMOLEC ) 

39 CONTINUE 

f 

C STEP-5- 

r COMPUTE initial VALUES OF LOCAL ULLAGE GAS VELOCITY 

PX a PCATAII) 

f. 
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V(N)«f LGM( 1I/AREA(N) 

NTEMP*N-1 

DO 40 L»1«NTEMP 

*i*K"L 

C10-Uh(J4n-»H( J} }*(2Uj4-l)42i(J)l)/(PX*(RA0(J-i-l}4-RA0( J> )*(CP( J4-1) 
1^CP( J)l) 

Cll s SURT(1*0>( (0R0X(J41) * 0R0X(J)>/2.Q)*«2) 

C12«0X*( ( OROXt J>1 )4-0R0XI ) / ( RA0( H-RA0( J) U 

Cl3*t »ZiU*l>**2+ZH J)**2)/JXJAY*PX*tZU'H» + Z< J))*(CPlJ+l»+CP( J) )l 
1 ) 

ci4-(Z2(j«n+z2( ji)/nz( j4-i)4-z (j) )*pxi 

C15-( (ZKJ+n-t-ZK J) )*(2.0«0IN>/( ( AREA { J4^1> 4-ARE A( J) )*PX*(CP ( J-4l}>CP 
l( Jl) ) I 

40 V( I ./< 1«-C12I*< ( 1«-4C12)*V< J*l)-C3*Cl0*{Tk( J+U4-T«( J)-Tt J+i) 

1 - T(J))4C11 -<2.*C3*C13 - DX*C14»40PDT - C3*C15) 

STEP-6- 

FINO GAS IN ULLAGE AT TIRE ZERO BY INTEGRATING DENSITY 


GASA=s(.5*AREA( 1)/ 1 T( 1 ) *Z ( 1 ) >-*- ( «5*AREA( N>) / I T (N )«Z ( N) ) ) 

NTEMP»N-l 

00 41 J»2*NTEHP 

GASA«6ASA«AREA( J)/T( J)/Z4 J) 

GASA»GASA«C4«P 

GSTART»6ASA 


WRI 


STEP-7- 

TE PROBLEM IDENTIFICATION AND INPUT DATA 


WRITE (6.1022) OT 


IF(UNUNS)250,251.250 
250 CONTINUE 

GSTART = GSTART ♦ •45354237 
WRITE (6.2015) GSTART 
GSTART * GSTART/. 45354237 
0^04 1054.3503 
WRITE (6.2050) C 
0 » 0/1054.3503 
TNKWT « TNKWT ♦ .45359237 
WRITE (6.4000) TNKWT 
60 TO 44 
2 51 CONTINUE 

WRITE (6.1015) GSTART 

WRITE (6.1050 ) 0 
WRITE (6.3000) TNKWT 

OIN - CINO(l) 

*0**m*4>^* *********** 
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r BEGIN MAIN PART OF CALCULATION 

C 

c STEP-a- 

C FIND TEMPERATURES AT NEM TIME 

C 

A4 CONTINUE 

CALL INTERP I TGAS .T IMEl. MTGAS. TI ME . TPI U ) 

CALL INTERP ( TSAT ,T IME^.MTSAT » TIME .TPI N*l» » 

CALL INTERP I T8ULK.T IME5. MTBULK.T I ME . TWPI N+l i I 
C 
C 

400 00 48 J-2.N 

D1=C5( J >*OT*OOUT/CPWI J ) 

02-l«0*C5(JI*DT4h( Jl/CPhl J» 

03=02/1 (C6/RA0IJ) )4H( JI*Z(J»*DT/P/CPI Jl» 

03 = 03/SORTIl.O^t(OROXl J>1)40R0XI 

04 s|C 7 *Z 1 ( JI^OPOT^-ICa/AREAlJ) M‘Z(JI*I-OIK)-*-C 8/XI N» *Z <J I / AREAIJ » 
1*(-0TR 1 >*OT/CPIJ ) /P 
r. 

60UA0=C3*I 1 .0+VCJ )*0T/0X-D4)-TM( J»-01 
COUAO=-C3*(TIJ»+VI J>*DT/0X*TPIJ-i> » 

RXl = -aOUAC/2. 

DISC = RX14RX1 - COUAO 
IFI0ISCI45. 45.46 

45 R0CT1= 0.5 ♦ ITU) ♦ TPU-D) 

TPU* = ROCTl 

GO TO 47 

46 RX2 = SORTIOISC) 

ROCTl = RXl ♦ RX2 

IF I ROOT D45C. 450. 460 
450 TPU) » 0.5 ♦ ITU) + TPU-D) 

GO TO 47 

460 TPU) = ROOTl 
C 

47 TWPU» = ITWU) + ID2-1.)*TPU)+D1)/D2 

48 CGNTINUE 

00 49 J=2*N 

01=C5U)«0T«0OUT/CPWU ) 

02= 1 .04C5I J )«OT«H U ) /CPw U ) 

49 CONTINUE 

C SPECIFIC HEAT OPTION FOLLOWING STATEMENTS 51.52 IS USED WHEN THE 

C TANK LIO IS (18-8) STAINLESS STEEL. THE EQUIVALENT THICKNESS 

C FOP THE LID MASS IS CONCENTRATED IN THE FIRST VOLUME ELEMENT. 

IFISPWSS - 499.0)53.50.50 

50 IFITWIl) - 75.0)125.51.51 

125 CPWII) = 0.010 

GC TO S3 

51 IFITwIl) .GT. 220. )G0 TO 52 

CPW(I)= 0.000416* TWIl) - 0.0203 
GO TO *3 

52 W = (TWIl) - 220.)/126.67 

CPtall) = 1(0.0018 ♦ W - 0.0127)*W ♦ 0.0374)*W ♦ 0.071 

53 C5( 1)=1.0/(TICK(1)«SPWSS) 

C9I 1)=XPI*0X«01A*T1CK(1)*SPWSS 
01=C5( 1 )*OT*OOUT/CPH( 1 ) 

02=1.0*C5( 1)*DT*H(1)/CPW(1) 
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c 


c 

54 

55 


C 

5l> 

C 

c 

57 

C 

c 

c 

c 

c 

c 

c 

c 


5 d 

152 

59 


60 

122 

C 


61 

151 

62 



T«P(1 )«4T>^m+(02-l«0J^TP( 1) ♦011/02 

CALL INTERP 4 OOUTO^T I ME6* MOOUT #T IMEtOOUT ) 
CALL INTERP 40IN0«TIHE7« MU1N*T1ME»01N1 
GO TO 54 

IF 4MPCATA-11 1«56«55 
PHCLO=P 

CALL INTERP 4 P0ATA.TIME2fHPCATA.TIMEtP 1 
TIPeP=TIME^DT 

CALL INTERP 4 POAT A«T IPE2 • HPOATA* Tl MEP« PPl 
£)PCT=4PP-P1/0T 

CCKTINUE 


IF 400001 73«73«57 
CONTINUE 

IF4 AB .EO* 0*1 GO TO 122 

STEP-9- 

EVALUATE ENERGY TRANSFER TO INTERNAL HARDWARE 


NOTE- CP-SPECF HEAT* 8K-PHEN0L1C# CU-COPPER* ST-STAINLESS* AL-ALUM 
PHENOLIC SPECIFIC HEAT DATA ESTIMATED FROM T PRC PUB* VOL 6 PT* II 
4 THERMOPHYSICAL PROPERTIES RESEARCH CENTE R-PUROUE UN1V*I 
00 60 J^ltN 

iF4TP4J) - TWB4 J)15£*5a*152 
T84J1 - TPlJl 
60 TO 60 

CP6K4J1 - 0*000664 ♦ TWBIJI 
TTT4 J1»TUB4 J1 

CALL HCOEFF 4 T4 J 1 *TTT 4 J) * P*XL *H4 J1 *PRAM4 J1 *ZCCNST*HMULT,HEXP»XMOLE 
ICl 

TB4J1*H4 J14^AB4t4T4Jl-TWe4 J) 1/CWTB«CPBK4 J1 l*0T^TliB4 J1 

001*001^H4 J1^AB*4TU1-TWB4 J11*0T 

CCKTINUE 

IF4A1C .EO* 0*1 GO TO 123 

CURVE FIT-CCPPER SPECIFIC HEAT DATA FROM WAOD TECH REPT 60-56 1960 
no 63 J*1*N 

IF4TP4J1 - TWC4 Jl>61*61*151 
TC4J1 = TPCJl 
GO TO Z3 

W ^ 4TWC4JI - 25*01/125* 

CPCU4J1 - 440*002I«W - 0*021^W ♦ 0*066831 ♦W 
TT14J1-TWC4 J1 

CALL HCOEFF 4 T4 J1 *TTT4 J1 *P *XL*H4 J) *PRAM4 J1 tZCONST*HMULT* HEXP*XMOLE 
1C) 

TC4 J1-HU1*A1C«4T4 J1-TWC4 J) ) / 4WT1C^CPCU4 J1 >^DT+TWC4 J) 

002=002^H 4 J 1«A IC* 4 T 4 J 1-TWC 4 J 1 1 *01 
CGKT INUE 

IF4A3S *E0« 0*1 GO TO 124 

CURVE FIT 418-dl STAINLESS STEEL SP HEAT FROM SCOTT CRYQ ENGR* - 
D* VAN* NOSTRANO 
00 66 J= 1 *N 
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IfUPlJ) - TSTIJ) >64*<4.135 
1^5 IF(TST(J» - 7i>.)130*131*131 

no CPST(J) s 0«01 

GC TO t5 

131 IFiTSTUl .GT. 220.01G0 TO 133 

f.PST(J) = 0*000«ia ♦ TSTIJI - 0.0203 
GO TO 65 

133 U > (TSTIJ) - 220.1/126.67 

CPST(J) * ( IO.OOia*M-0.0127)6M ♦ 0.0374)*W ♦ 0.071 

64 T71J1 3 TPIJI 
GC TO 66 

65 TTT(JI*TSTU) 

CALL FCOEFF 1T( J 1 .TTT I J 1 . P .XL .HUl .PRAM ( J I .ZCQNST. HMULT, HEXP . XMOLE 
If.) 

T7( J)«H(J)*A3S*(TU)-TSTIJ) )/( MT3S6CPST( J) )*DT+TSTIJ) 

003s003+h< J l*A3S*m J l-TST ( J 1 1 *DT 

66 CChTINUE 

124 IFlASe «E0. 0.1 GO TO 126 

C CURVE FIT GF AL . ALLOY 6061-T6 SPECIFIC HEAT DATA FROM TPRC 

OC 72 J*1.K 

IFITPIJl • TALStJ))70*70. 153 

70 TBCIJl s TPfJ) 

GC TO 72 

153 IF(TALSIJ) > 70. 0167.6a. 68 

67 CPAL(J) - 0.000397 * TALSIJl - 0.013 
IF(TALSIJ) .LT. 36.01CPALI Jl-0.0012 
GO TO 69 

68 H - ITALSIJ) • 70.1/117.5 

CPALIJ) = I(0.00334*M • 0.03511«M * 0.13666)«M * 0.015 
IF(TALSIJ) .GT. 540.) CPAL(J) * 0.212 

69 CCNTI6LE 

71 TTTI J IsTALSIJ) 

CALL HCOEFF I T( J 1 *TTT I J1 . P .XL. HlJ 1 « PRAM! J) . ZCONST . HMULT. HE XP. XMOLE 
ICl 

TBCIJl’^HI J)*A5e«(T( J)>TALS( J)1/(WT5B6CPAL(3) 1*0T>TALS( J) 

004=OU4«H( J1«A5B*(T( J1*TALS( J) 1*0T 

72 CCNTINUE 

126 OIN = 4001 ♦ 002 + 003 * 0041/ 4T1ME * X(N1) 

73 CONTINUE 
C 

r STEP- 10- 

C FIND COMPRESSIBILITY FACTORS AT NEW TIME 

C 

IF (2CCNST) 1.74.76 

74 ZHCL0=ZI51 
00 200 K=1.N 

200 ZAVIKl = 2IK) 

OC 75 J=1.NP 

75 CALL CCHPRS I TP IJ) . P .Z O 1 . 21 ( J 1 . Z2 I J) . XMOLEC 1 

76 CONTINUE 
( 

C 
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c STEP-ll- 

C fINO VELOCITIES AT NEIi TIME 

C 

VHOLOsV<5) 

CALL INTERP ( FLCkt* T IME3, MFLOW , TI ME*FLQWNP > 
V(NP) s FLCMNP/AREAINPI 


C 

77 


C 

C 

C 

C 

C 


78 

C 

c 

c 


79 

C 

C 

C 


C 


C 

c 

c 

80 

c 

c 

c 

c 

c 

61 


300 


6P 

63 

84 


no 77 J=1*N 
K*hP-J 

V(K) = (TP(K)*V(K+l)-(OX/(OT*Z(Kn»*(Zl(KI*(TP(K)-T(K) n 
I -f Z 2 <K>*TP(KI*DPDT#DX/I Z(Kl*PPn/lTP(K»+(illK)/Z(K)>*(TPIK+i) 
? - TPIK)» - I2.0*OX ♦TPIKl ♦ IOKOX(K+11+OROX(K>)»/(RAO(K+U 
3 * KACIK))) 


STEP-12- 

FINC HEAT FLOW RATE TO MALL AND TOTAL HEAT AOOEO TO HALL 

00«C9( 1 )*CPM(1IA(TWPU)-TM( in 
00 78 Js2*N 

00>004RA0( J)*S0RT(1«0+0R0X(JH<*2)*CPM( J)«( TMPIJ1>TM( J) )AC9(J1 
UsC4D0 


FIND GAS FLOW RATE AND TOTAL GAS AOOEO UP TO THE NEW TIME 
GASBs0.5*AREA( U/TP( 1)/Z( mo. S*AREA(NP) /TPI NP )/Z ( NP> 

00 79 J«2.N 

GAS8-GAS64AREA(J}/TPU1/Z(J> 

GASBaGAS6*C4*P 

NOTE GASCHECK CALCULATION BASED ON CROSS SEC ION AT NET POINT »5 

GSRATEsABSIIGASB-GASAl/OT) 

GAS^GASB-GSTART 

GASCHK*GASCHK40.25*I VI 51 ♦VHOLO )/ (C8/AREA (5)1 *DT*( PHOLO/ZHOLO/T (51+ 
IP/Z(51/TPI511 


STEP-13- 

FINO SPECIFIC HEATS AT NEW TIME 
00 80 J»1*NP 

CALL SPHEAT ( TP( J 1 .TWP U 1 1 CP ( J 1 tCPWlJ 1 . XMOLEC) 


STEP-14- 

FINO HEAT TRANSFER COEFFICIENT AT NEW TIME 

IF(TACC182.82.81 
XL»2.0«RAO(N41 1 
TTT(N4l)aTP(N+l) 

CALL HCOEFF ( T AOO.TTT ( N4 1 1 .PtXLtHiNi-l) ,PRAN( N^-ll .ZCONST.HMULT t HEXP 
I.XMOLECl 
H2H«H(N4ll 

HSUR = H2H * .14/HMULT 

OTR<HSUR«AREA( N4l )«( TAOO-TP (N4 1) 1 

C0TR«C0TR40TRW0T 

GC TO 83 

UTR=0.0 

IF (HCCNSTl 1.84*89 
CONTINUE 
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nr as j-i*np 

If (X( JI-RCONSTI 85«£5«86 

C IN TURBULENT RANGE HEXP - 0«333t CHOICE OF XL IS IMMATERIAL 

£6 Xl^CIA 

GO TC €7 
£6 XLsCIA 

£7 TTTIJ»aTRP(J» 

CALL hCOEFF ITP <J » .TTT ( J I , P .XL .HIJ I , PRAM! J > . ZCONST ,HMULT t HEXP, XMOL 
lEC) 

IF (TR2.EC.O.) GC TO 88 

WBEXPO = EXPl-0.00117 ♦ CARAD**2*X (Jl ) 

WhSf!*01C/CA8AD**.2*IGSRATE/ADIFUI**.8*PRAM( J I 

88 HI J)sH(J)*TRl'*'8HSO«MBEXPO*TR2 
C 

89 CONTINUE 
C 

C 

ITERl * ITERl ♦ 1 

IF I ITERl'lWRiTE) 92,91.91 

91 CALL NRITE2 
ITE«1*0 

C 

C 

C STEP-15- 

C CHECK TO SEE IF END TIME HAS BEEN REACHED 

C 

92 IF ITINE-ENDTIM) 93,98,98 

93 TVME=TIME*OT 
C 

C END TIME NOT REACHED - PREPARE FOR ANOTHER STEP 

C 

CALL INTERP ( FLOh ,T IME3, MFL08, TTME ,FL08P) 

VPsFLONP/AREAINP) 

OT » DX/IWINP)4^VP> * 2, 

TIME=TIME*OT 

C 

DO 9A J-l.NP 
Tl JI»TFIJ> 

94 TtalJI^TWPlJl 
IFICtiOC>97,97,95 

95 DO 96 J»1,N 
TUBI JlsTBIJI 
TmCIJ»=TCI J> 

TSTIJlsTZIJl 

96 TALSI J)-TBCIJ) 

TtoB(N*ll-T8PIN4-U 

TWCIN411*T8PIN4^1I 

TST(N«1)sTWP(N4-1) 

TALSIN-»l>-TWP(N«ll 

97 CONTINUE 
f 

N-N4^1 

NP«NP^1 

CNsCN-i-1.0 

GASA^GASB 

C 


GC TO 44 


40 


c 

c 

c 

c 

c 

c 

9H 

C 


99 

C 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

100 

101 

102 

103 

1015 

1022 

1050 

2015 

2050 

3000 

4000 


STEP-16- 

END TIME EXCEEUEC - INTERPOLATE CONDITIONS AT END TIME 


RATIO-CENCTIM-TIME^OT>/OT 

TIME^ENOTIM 

no 99 J«1«N 

TPCJ1 = TIJI+RATI04(TPU1-T( J1 1 
T^PI J)-TU( JI>RATI04ITUP(JI-Ti«( J) I 

O-C-OC^RAT 10400 

GAS»GAS-GASB-fGASA4RATI04IGAS6-6ASA) 
GASCHKsGASCHK-GASB^GASA^RATI04IGAS&-GASAI 
XI NPMXCN)4RATI040X 

STEP-17- 

liRITE OUT RESULTS REFER TO SUBROUTINE MRITE 2 
CALL NRITE2 

GO TO 1 




FORMAT STATEMENTS 


FORMAT 

FORMAT 

1 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

END 


(lHlt30X«24H TANK EXPULSION PROGRAM/ IHJI 
I80H 

) 

(8F10.0I 

(FlO.Ol 

(lHKf23H INITIAL ULLAGE GAS ^ F6«3f5H LBS) 

(1HK« 27H INITIAL TIME INCREMENT « F6«It9H SECONDS) 
UHK^ 25h INITIAL HEAT TO mALL = F7*lf 5H BTU) 
(1HK«23H INITIAL ULLAGE GAS » F6«3«10H KILOGRAM) 


UHK* 

25H 

INITIAL 

HEAT TO 

taALL 

= Ell 

• 4* 7H 

JOULE) 

( IHK, 

29H 

THE 

TANK 

MEIGHT 

LESS 

LIO a 

F7.1, 

5H 

LBS) 

( IHK « 

29H 

THE 

TANK 

mEIGHT 

LESS 

LID s 

F7.<i, 

lOH 

KILOGRAM) 
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TANK EXPULSION PROGRAM 
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APPENDIX C 


LISTING OF RAMP PROGRAM 


c 

c THIS IS THE PRESSURIZATION MAIN PROGRAM 

♦♦♦♦♦♦♦ 

DIMENSION TGAS(IOO), TIMEl(lOO), POATAdOOl , TIMt2(100), FL0WI25I, 
1 TIME3{25), TSAT(50», TIME4(50I, TBULKI25) , TIMES 125), COUTDI25), 
2TIVIE6I25), 0IND(25), TIME7(25), ZAVI250) , TWINDIIOO), DISTKIOO), 
3TGIN0(1C0), DIST2UC0), XI250), TI250) , TPI250), TWI250), TWP(250) 
4, V(25C), VP( 250), CP(250), CPWf250), H(250), ZI250), ZK250), Z2I 
5250), RAD(250), AREA! 250) , YRA0(250), XRADI250), DRCXI250), TtST(2 
650) 

C 

DIMENSION C5(250), C9(250), YTICKI250), XXTI250), TICKI250), XR(3) 
1, A0(3) 

COMMON X,T,TP, TWfTMP, V,CP ,CPW,TGAS, 71 Me I , MTGAS , POAT A, T I ME2 , MPDAT A, 
IFLQW, TIME3,MFLOW, TSAT , TI ME4 ,MTSAT, TBULK , T I ME5 , MT8U LK, OOUTO , TI ME6, M 
200UT,Q INO, TIMET, MQIN,N,NP ,XN .ULLAGE , RAD IUS,TL ID ,SPVIGT ,XMOLEC, GSTAR 
3T ,HC0NST,TIME,GAS,O,GASCHK,GSRATE ,0T, ZCONST ,H ,HMULT ,HEXP,YRAO,X RAD 
4, MRA0,DQ,0IN,UNUMS,TWIND,DISTI,MTWI N,TGINO,DI ST2, MTGI N ,DIA, 

5TA0C, YT ICK,XXT,MTICK 

♦♦ 


READ 3 CAROS OF PROBLEM DESCRIPTION AND WRITE OUT. THERE 
MUST BE THREE CAROS USED, ALTHOUGH ANY OR ALL OF THEM MAY 
BE BLANK. LEAVE THE FIRST COLUMN OF EACH CARD BLANK AND 
ENTER ANY INFORMATION IN COLUMNS 2 TO 30. 

WRITE (6,106) 

DO 2 J = l,3 
READ ( 5, 107) 

WRITE (6,107) 

This IS THE PRESSURIZATION PROGRAM 


NOMENCLATURE FOR INPUT DATA 


XN, NUMBER OF NETPOINTS AT TIME ZERO (MUST BE 3.0 OR MORE) 

OUTPUT, NUMBER OF TIME STEPS TAKEN BEFORE EACH OUTPUT 

ULLAGE, INITIAL ULLAGE LENGTH, FEET (CANNOT BE ZERO) 

RADIUS, TANK RADIUS, - SPECIFY WHEN TANK IS A CYLINDER 
THE PROGRAM THEN ASSUMES SPHERICAL END SECTIONS. 

IF RADIUS = 0., THE RADIUS IS INTERPOLATED BASED CN THE 
TANK RADIUS (VS DISTANCE FROM TOP) DATA READ INTO THE 
PROGRAM . IF RADIUS = -1, THE TANK IS A SPHERE. 

SPWGT, TANK WALL SPECIFIC WEIGHT 
SPWSS, SPECIFIC WEIGHT OF TANK LID MATERIAL 
ENDTIM, TIME TO COMPLETE THE HOLD PERIOD 
XMOLEC, MOLECULAR WEIGHT OF PRESSURIZING GAS 

ZCONST, COMPRESSIBILITY FACTOR (1. FOR IDEAL GAS, BLANK FOR REAL) 
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c 


HCONST, HEAT TRANSFER COcFF. IBLANK IF H IS TO BE CALCIJLATEDI 
BTU/HR, /SQ.FT. /DEG. F. ♦tWATTS/SO. M/DEG K 
DT, INITIAL TIME STEP (SECONDS) 

ENDRMPf TIME TO COMPLETE THE RAMP PERIOD (SECONDS) 

OIA» DIAMETER WHEN THE TANK IS A SPHEREt OR CYLINDER. 

SPECIFY DIAMETER ON VERTICAL AXIS WHEN TANK IS SPhERCID 
RlfP2t TTME COORDINATES! SECONDS) ALONG THE PRESSURE RISE CURVE 

P3,P4, THESE CONTROL THE INITIAL SELECTION OF THE TIME STEP 

SEE EFN STATEMENT 96 ET CETERA FOR SIGNIFICANCE 
CYLN, WHEN RADIUS IS SPECIFIED* THE CYLINDRICAL LENGTH BETWEEN 

THE TANK END SECTIONS MUST BE DEFINED. 


IF HCONST IS BLANK H WILL BE COMPUTED FROM THE EQUATION 
H = HMULT*CONO/XL*(NUSSELT*PRANDTL)**HEXP 

HMULT, CONSTANT IN ABOVE EQUATION (0.13 IS USED IF LEFT BLANK) 

HEXP, CONSTANT IN ABOVE EQUATION (0.333 IS USED IF LEFT BLANK) 

TLIO THE MASS OF THE TANK LID AND FLANGE MADE EQUIVALENT 

TO THE FIRST WALL THICKNESS 

UNUMS* SPECIFY GREATER THAN 0. TO PROGRAM IN S I UNITS 


CROUP I DATA RAMP DATA 

*ltirHLi^*THi4t*Jtt^i*4iitl*:4fitLm************ 

NAMaiST / PARAMS / XN .OUTPUT * ULLAGE^RAOI US *S PWGT * E NOT I M,X MOL EC, 
1HCONST,DIA,ZCONST,HMULT,HEXP,SPWSS,OT,ENORMP, PMl,PM2, PM3,PM4,CYLN, 
2UNUMS 

READ (5, PARAMS) 


CROUP II DATA RAMP DATA CONTINUED 

#]|c # 3|ci|u|t I|( Itc # 


NAMELIST / INPNTS / MTGAS.MPDATA.MFLOW.MTSAT, MT8ULK,MOOUT, MQIN. 
IM TW IN , MTG IN , MR AD , MT IC K 
READ (5, INPNTS) 


NAMELIST / OTPNTS / MTGAS.MPDATA.MFLOW.MTSAT, MTBULK.MQOUT, MQIN, 
IMTW IN.MTGIN.MRAO.MTICK 
WRITE (6, OTPNTS) 

NAMELIST / TABLES / 


1 TGAS , TIMEl . PDATA 

2 TSAT , TIME4 . TBULK 

3 QINO , TIMET , TWIND 

4 YRAD , XRAO , YTICK 

READ (5, TABLES) 


T1ME2 
TIMES 
01 STl 
XXT 


FLOW 
QOUTO 
TGI NO 


TIMES 
TI ME6 
0IST2 


PRELIMINARY COMPUTATION 

GSTART=0.0 
T IMY*0.0 
LOOM*0 
ITERI=0 
0 * 0 . 

GASCHK*C. 

TICK! 1)*YTICK( 1) 

TLIO*TICK( 1) 

CALL WRITE 1 
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c 

301 

r, 

c 


3 

4 

5 

6 

7 

8 

9 


IFfiJNUMS .EO. 0. I GO TO 301 
CALL SINTS 
CONTINUE 
TPP=l . 

XPI=3 .14159 
R=1545 .4 
X J4Y=778.2 
HC0NST=hCONST/3600. 


N =XN 
NP=N + 1 

DX = ULLAGE/(XN-1.0) 

IWR ITE=OUTPUT 
CHECK = FLOW( 1) 

0ETY=07 

C8=R/XM0LEC 
C3=CX*C8 
C l=C3/2. 

C2=C1/XJAY 

C4=CX/C8 

C6=2.*C8 

C7=C8/XJAY 

UWTP = 2.0*XP I*DX*SPWGT 
on 3 J=2t250 
X JTICK=J“1 
X(J »=DX*XJTICK 

CALL INTERP ( YTICK , XXT ,MTIC K, X( J ) 
C5( J)=1.0/(TICK( J )*SPWGTI 
C9(J) = UWTP * TICKU) 

X(11 = C. 


STEP- lA- 

MAKE GEOMETRY CALCULATIONS 

IF (RACIUSJ 12f4,6 
DO 5 J=2,250 
XJTEMP=J-1 

CALL INTERP ( YRAO ,XRAD ,MR AD , X ( J > 
AREA! J ) = XPI*RAO(J )♦♦2 
R ADI 1 l=PAni 2) 

AREA! 1)=APEA( 2) 

GO TO It 

CONTINUE 

DO 10 J=2,250 

IF (XI JJ-DIA/2.) 7,8,8 

RADIO >=SORTIOIA#XIJI-Xt 

AREAI J )=XP I*RADI J »**2 

GO TO 10 

IF IXI J).GT.ICYLN+RAOIUS» > GO TO 

RADIO |=RADIUS 

AREAI 0 ) = XPI*RAOIO )**2 

GO TO 10 

PL6G=XI0 l-CYLN- RADIUS 


,TICK(OI ) 


RADIO) ) 
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10 

11 


c 

12 


13 

lA 

15 

C 

16 
C 
C 
C 

C 

C 

C 

C 

c 


c 

c 

c 

c 

c 

17 

c 

c 

c 

c 

c 

c 

18 

19 

C 

20 
21 
C 

22 

23 
C 

24 

25 

26 


IF (PLEG.GE.RADIUSI GO TO 11 
RAOU ) = SQRT(RADIUS**2-PLEG**2> 

ARF A( J )*XPI*RAD(J »**2 
CONTINUE 
R ADC 1 l*RAD( 21 
AREAC 1 )=AREA( 21 
GO TD 16 

00 14 J*2f250 
X JTEMP=J-l 

IF (X( J J-DIAl 13,13,15 
RADCJ > = SQRTC0IA*X(J1-XU1**2I 
AREAC J l«XPI*X CJ )*<OIA-XC Jl I 
AREAC 1)*AREAC 21 
RAOCl )=RA0C2) 

CONTINUE 

THE NEXT 3 CAROS ARE SKIPPED SINCE THE FLOW DATA IS ZERC IN RAMP 
DT=DX*AREA(N l/FLOWC 1) 

CALL INTERPCFLOW,TIME3,MFLOW,DT,FLO) 

DT* DX* AREACN)/CFLOW(1H-FLO)*2.0 
TIME=CT 


STEP-IB- 

COMPUTE INITIAL WALL AND GAS TEMPERATURES 
00 17 J»1,N 

CALL INTERP C TWINO,DI ST1,MTWIN,XIJ) ,TWC Jl) 

CALL INTERP C TGINO, 01 ST2, MTGIN, XC J) , TC Jl ) 


STEP -IC- 

COMPUTE INITIAL VALUES (DF SPECIFIC HEAT 
CALL SPFEAT C TC J 1 ,TWC J 1 ,CP C J) ,C PWC Jl , XMOLEC I 


STEP- 2- 

COMPUTE INITIAL VALUES OF HEAT TRANSFER CCEFFICIEKT 

IF CHCONSTl 1,20,18 

DO 19 J»l,250 
HCJ )=HCONST 
GO T3 26 

IF CHMLLTl 1,21,22 
HMULT=C.13 

IF CHEXPl 1,23,24 
HEXP=0.333 

DO 25 J*l,N 
XL=ULLAGE+XC J 1 

CALL HCOEFF C TC J1,TWC J) ,POATA C 1 1 ,XL,H( J 1 , ZCONST ,HMULT , HEX P, XMOLEC) 
CONTINUE 
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r 


c 

c 

c 

c 


27 


28 


C 

29 

30 

31 
C 

c 

c 

c 

c 


c 

32 

c 

33 

34 

C 

35 


C 

c 

c 

c 

c 

36 


STEP- 3- 

CQMPUTE INITIAL VALUES OF COMPRESSIBILITY FACTOR AND 

ITS DERIVATIVES 

IF (ZCCNST) 1,29,27 

no 28 J=l,250 
7 (J )=1.0 
Z 1( J)=l .0 
Z2( J)=1.0 
ZHOLO=1.C 
GO TO 21 

DO 30 J = 1,N 

CALL COMPRS ( T( J » ,POA TA ( 1 1 , Z ( J) ,Zl ( Jl ,Z2 (J» ,X MOLECI 
CONTINUE 


STEP-4- 

INITIALIZE VALUES FOR OOUT, QIN, P AND OPDT 

00iJT=Q0UTD( 1) 
g IN =-OINDI 1) 

P=P DATA! 1 > 

PKDLD=P 

IF (MPCATA-U 1,32,33 
PP=P 

GO TO 24 

CALL INTERP ( PDATA,TI ME2, MPOATA , TI ME , PP) 

OPOT=( PP-P 1/OT 
OPDTPV*OPDT 

00 35 J=2,N 

ORDXI J l*(RAD( J*ll-RAO( J»l / (XIJ<-l)-X( J1 1 
0R0X(N*1 )=DROX(N) 

ORDXI I )=(RAO( 2I-RADI 11 1/XI2I 

STEP- 5- 

COMPUTE INITIAL VALUES OF LOCAL ULLAGE GAS VELOCITY 

VIM)=FLCWI ll/AREAIN) 

NTEMP=N-1 

DO 37 L=1,NTEPP 

J=N-L 

C12=DX*( (DRDXIJ+lKOROXUn /(RAO(J+ll+RAD( jn I 
C10=H(J )*Z1(JI/(P*RA0(J)*CP( J)l 
C11=SQRTI 1.0+(I0RDX(J«-ll+DR0X(jn/2.»**2l 

VA=C3*C10*(TWIJ l-TIJl l*Cll+ ICl /AREA! J)) ♦Zl ( J» ♦CIN/ P/CP ( Jl* I C2*Z lU 
l)**2/CP(J »-.5*DX*Z2(J) l*DPDT/P/Z(J) 

VA=VA/( 1.0-C12) 

K=J*l 
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37 


C 

c 

c 

c 

c 


38 


C 

C 

C 

C 

c 

c 

c 

c 

250 


251 

C 

C 

120 

C 

c 

c 

c 

c 


c 

c 

c 

c 


39 

c 


C20*H( K )*Z1(K )/(P*RA0 (KI*CP(K>) 

VB=C3*C2C*( TW(K )~T( K))*C !!♦ (C 1 /AREA (Kl ) *Zl ( K) *QIN/ P/CP ( KH- ( C2*Z 1( K 
1) **2/CP(K >-.5*DX*Z2(K) )*DPDT/P/2(K) 

VB=VB/( 1.0-C 12) 

V(J ) = ( 1.0+C 12)/U.0~C 12)*V( J+1)-VA“V8 
WRITS (6,111) 

WRITS (6,112) ( V( J),J=1,N,1C) 

IF (GSTART.GT.O.) GO TO 40 

STEP- 6- 

FIND GAS IN ULLAGE AT TINE ZERO BY INTEGRATING DENSITY 


GASA=0.5*AREA( 1)/T( 1) /Z ( 1 ) +0. 5*AREA( N) /T ( N) /Z ( N ) 
NTEMP =N-1 
00 38 J=2,NTLMP 
GASA=GASA4APEA( J)/T( J)/Z( J) 

GASA=GASA*C4*P 

GSTART*GASA 


STEP- 7- 

WRITE PROBLEM IDENTIFICATION AND INPUT DATA 


I F( UNUMS )250, 251, 250 
CONTINLE 

GSTART = GSTART * .4535924 
WRITE (6,2015) GSTART 
GSTART = GSTART/. 4535924 
Q = Q 4 1054.3503 
WRITE (6,2050) Q 
0 = 0/1054.3503 
GO TO 120 
CONTINLE 

WRITE (6,1015) GSTART 
WRITE (6,1050 ) Q 
CONTINUE 

*41 * ^ 4c 4t 4 4c 4c » 4c 4c 4c 4c 4c 4c 4c 4e 4t 


EEGIN MAIN PART OF CALCULATION 
T( 1 )=TGAS( 1) 

INITIALIZE AN ESTIMATE FOR TP(J) 

NJ=N- 1 
00 39 KI*1,N 
TP(KI )=T(KI ) 
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STfP-8- 

FIND TEMPERATURES AT NEW TIME 


C 

C 

c 

40 


C 


C 

r 


41 


c 

42 

C 

C 


43 

44 

45 

C 

c 

46 
C 


C 


C 


C 


CALL INTEPP ( TGAS,TIMEl,MTGAS,TIME,TP(l) I 
CALL INTEPP (TSAT,TIM£4,MTSAT,TIMEtTP(N»» 
CALL INTEPP ( TBULK,TI ME5,MTBULK,TIMfc,TWP(N») 


K ST AR = C 
7 HnL0 = 2< 2 ) 

VH0Ln = V( 2 ) 

FIND MISCELLANEOUS QUANTITIES AT NEW TIME 


CALL INTEPP 
CALL INTERP 
OIN = -OIN 
IF (MPCATA-l) 

PHOLO = P 
CALL INTERP 
TIMEP=TIME+DT 

CALL INTEPP (P0ATA,TIME2,MPDATA,TIMEP,PP) 
DPOT=(P-PHOLD)/DT 


(00UT0,TIME6,MQ0UTt TIME tQOUT) 
(QIND,TIME7,M0IN,TI MEfQINI 


l»42,41 


(P0ATA,TIME2tMP0ATA,TIME,P» 


CONTINUE 


INITIALIZE AN ESTIMATE FOR VELOCITY 
DO 45 J=2»NJ 
I F ( DPCTPV-0. ) 43f 43, 44 
V(J ) = V( J ) 

GO TO 45 

V(J I*V(J »*DPDT/DPDTPV 

CONTINUE 

DPOTPV=DPDT 

KLAMP=0 

I STAR =C 
CONTINUE 
00 54 J*2,NJ 

D1=C5( J J*DT*QOUT/CPW( J» 

02=l.0*C5( J )*DT*H(J)/CPW( J) 

D3=D2/< (C6/RA0( J) »*H( J )*Z U>*DT/P/CP( J)) 

D3=D3/S0RT( 1.0+ ( ( ORDX ( J+1 H-DRDX( J) » /2.» ♦♦2) 

04={ C7*Z1( J )*DP0T+(C8/AR£A( J) )*Z( J) ♦QINI *DT/C P( J)/P 
G 1=D3-C3*D4-TW( J 1-0 1 
G2=C3*CT/DX 
G3=D3*T( J ) 

G4=71( J )/Z( J»*DX/DT 
G5=Z2( J»/Z( J )»DX*DPDT/P 

G6=G5-C4 
G7=G4>CT( J I 

G8=1.0-Z 1( J >/Z( J)-2.*0X/RAD ( J)*dROX( J l 
PHE=Z 1( J )/Z(J )*TP( J+1H-G1*G 8*G2*( V( J^-1) *G6) 


52 


IIU=GI*ZI(J)/Z ( J f*TP( J*1)<-G2*G7-G2*TP( J-1)*(V( J+ 1 »+G6) -G3*G8 
AAY = l./2.»( 3.*U(J-PHE**2» 

RR=-G2*G7*TP( J- 1 J-G3* Z 1( J ) /Z ( J) *TP( 1 1 
8Ee = l*/27.X'(2.*PHE**3-9.*PHE*UU+27.*RR» 

C 

0 IS l=Ut**2-4.*PHE*RR 

IF (DISI.LE.O.I GO TO 47 

TP( J)=(-UU*S0R7<DIS1» l/(2.*PHE> 

GO TO 48 

47 TP( J»=C.5*(T( JI*TP(J-1I) 

48 DO 49 LEL=1, 10 

HOP E= G 8*TP( J )♦♦ 34-PHE* TP ( J »**2«-UU*TP ( J) +RR 
HOP =ABS( HOPEI 

0K)PDT = 3 .*G8*TP(J l**2+2.*PHE*TP< JUUU 
AZ9=800. 

AZ8 = .0C01*0H0PDT 
IF <AZ8.LT.AZ9I AZ9=AZ8 
IF (H0P.LT.AZ9) GO TO 50 
TP( J|*TP( JI-HOPE/OHOPOT 
XR( 1I=TP( J » 

49 CONTINUE 
WRITE (6,105) 

KSTAR=KSTAR*1 

IF (KSTAP .GT.25) GO TO 55 

50 CONTINUE 

TEST( J )=BEE**2/4.+AAY**3/27. 

IF (TEST(J).LE.O..ANO.TP( JI.GT.TPP) GO TO 51 
GO TO 54 

51 A0(l)*C-8 

A0( 2)=PHE*XR( 1)*G8 
A0( 3)*UU*XR( 1)*A0(2) 

XQl=-AQ(2)/( 2.*AQ(in 
OISC=XCl*XO 1-AQ(3)/AQ( 1) 

52 IF (DISC .LT. 0.) GO TO 54 
C 

53 X02=S0RT(DISC) 

TP( J)*X01+XQ2 

IF (TP(J).GT.TPP) TP(JI=XR(1) 

C 

C FOR SIGNIFICANCE OF TEST(J), SEE COMMENT UNDER STEP 13 

54 TWP(J »=(TW(J)*(02-1.0)*TP(J) + 01)/02 
GO TO 56 

55 IF (TIME.GT.PMl) DETY*OETY-0. 005 
IF (TIME.LE.PMI) OETY=0ETY-0. 01 
DT=D£TY 

TTME=TIMY+OT 

IF (TIMfc.LE.TIMY) GO TO 1 
GO TO 36 

C THE FOLLOWING TWO EQUATIONS ARE SPECIFIC HEAT VS TEMP. FOR AN 

C ASSUMED MASSIVE STAINLESS STEEL LIO AT THE TOP CF THE TANK 

56 IF (SPVSS-499,0) 61,57,57 

57 IF (TW(l)-75.0) 58,59,59 

58 CPW(l)=C.C10 
GO TO 61 

59 IF (TW( 1 ).GT.220. ) GO TO 6C 

CPW( I ) = C.0004ia*TW( 11-0,0203 

GO TO 61 
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60 W=( TW< n-220. )/l26.t7 

CPW( n=( (0.00 18*W-0. 0127)* W+0. 0374) ♦W+0. 071 

61 C5( 1)=1.C/(TICKU)*SPWSS> 

C9{ l)=XP I*0X*0IA*TICK( U*SPWSS 

c 

Dl=C5( 1 )*DT*00UT/CPW( 1) 

0 2=l.0-»C5( l)*DT*H( l)/CPW( 1) 

KLAMP=KLAMP+1 

c 

TWP( I ) =(TW( I )+(D2-1.0)*TP( l) + Dl) /02 

C 

IF (KLAMP.lt. 3) 60 TO 62 
IF (ERPP .LT..05) GO TO 62 
IF (ERPP.GT..05.AN0.ISTAR.GT.6) GO TO 62 
GO TO 71 
C 
r, 

C STEP-9- 

C FIND HEAT TRANSFER COEFFICIENT AT NEW TIME 

C 

62 I F (HCONST) l,63f 65 

63 DO 64 J=l,N 
XL=ULLAGE*X( J ) 

64 CALL FCOEFF ( TP (J ) , TWP ( J ) , P , XL» H ( J ) , ZCONST ,HMULT *H EXP ,X MOL EC > 

65 CONTINUE 
C 

c 

C STEP- 10- 

r FIND COMPRESSIBILITY FACTORS AT NEW TIME 

C 

00 66 K=1,N 

66 ZAV(K)=Z(K) 

I F (ZCCNST) 1,67, 69 

67 00 63 J*1,N 

68 CALL CCMPRS ( TP ( J ) ,P, Z ( J) , Z1 ( J) ,Z2 ( J) , XMCLEC) 

69 CONTINUE 
C 

C 

C STEP- 11- 

C FIND SPECIFIC HEATS AT NEW TIME 

C 

00 73 J=1,N 

70 CALL SPEEAT ( TP ( J ) , TWP ( J) ,C P ( J) ,C PW ( J ) , XMOLEC ) 

C 

c 

71 CONTINUE 

1 F (KLAMP.lt. 3) GO TO 46 

C 

C STEP- 12- 

C FIND VELOCITIES AT NEW TIME 

C 

CALL INTERP ( FLOW ,T IME 3 ,MFLOW , T I ME , FLOWNP ) 

V(N ) = FLCWNP/AREA(N) 

I F (LrCM .GE .3 ) TO 89 
C 
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72 


73 


74 

75 

76 


77 

78 

79 
83 

81 


82 

83 

84 

85 

86 


C 

C 

87 

88 

C 

C 

C 

c 

89 

C 


on 72 J=2tN 
K =MP-J 

V(K ) = ( TP (K>*V(K + U- (DX/(OT*Z(K) ) )*(Z1 (K)*(TP( K)-T(K) ) H-Z2( K)*TP(K) 
1#DPDT*CX/(Z{K )*PP ) )/( TP(K» + (Zl( K» /Z(K))*(TP(K*-1 )~TP<K> ( 2 . O+QX *T P 

2 ( K )*( CRCX(K+l )*DRDX(K n ) / (RAO(K*U4-RAD(Kn ) 

I F ( V (K ) .LT.O. > GO TO 73 

COMTINLE 

ERRP=C. 

GO TO 83 
LOnM=LCOM,+ l 
WRITE (6,119> 

GO TO 55 
DO 75 J=l»N 
VP( J )= .5 
CONTINUE 
CONTINUE 
DO 78 J=1,N 
VIP=VP(J ) 

IF (V IP .to. 0. 1 GO TO 78 
ERR =( V( J )-VP( Jl )/VP (J ) 

ERR=ABS( ERR ) 

I F ( ERR-ERRP ) 78, 78,77 

ERRP=ERR 

CONTINUE 

IF (ERRP-.004) 89,89, 79 

IF ( ISTAR ,E0,36) GO TO 87 

DO 81 J = 1,N 

VP( J)=V( J I 

CONTINUE 

I STARS ISTAR + l 

I F (ISTAR .GT.32 ) GO TO 103 
IF { ISTAP-40) 46, 46,84 
I F ( I STAR ) 74, 74, 76 
IF (ERRP-.IOI 85, 85,86 
ERROR =ERRP*100, 

GO TO 10 4 
WRITE (6,1081 
WRITE (6,114) ERROR 
WRITE (6,115) TIME 
WRITE (6,110) 

WRITE (6,112) (TP(J), J=1,N,10) 

write (6,116) 

WRITE (6,112) (TWP( J) ,J=l,N,10) 

AT THIS POINT THE PROGRAM WILL CONTINUE REGARDLESS CF THE ERROR 

OPTITN II IS TO REDUCE THE TIME INCREMENT, PROCEED TO STEP 5 

GO TO 89 

DO 88 J=1,N 

V(J )=( V( J )+VP( Jl) /2. 

GO TO 8C 

STEP-13- 

FIND GAS FLOW RATE AND TOTAL GAS ADDED UP TO THE NEW TIME 
CONTINUE 

IF (TP( D.GT.TPP) TPPsTP(l) 

IF TEST(J) IS GREATER THAN O-ONE REAL AND TWO CONJUGATE 
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c 

c 


90 

91 


92 

C 


c 

c 

c 

c 

c 


93 

c 

c 

r 

C 


94 

95 

C 

C 

c 

c 

c 

c 

c 

c 

c 

96 


lyAGINERY ROOTS — IF T£ST(J)=0. t THERE WILL BE THREE 
REAL ROOTS — IF TEST(J) IS LESS THAN D-THGRE WILL 8E 
IFPLE REAL AND UNEQUAL ROOTS, 
no 91 j=i,N 

I F (TPU HTPP ) 91,91t 90 

WRITE (6,109) TP( J), JfTIME ,TEST(J) 

CONTINUE 

GASB=C .5*AREA ( 1 ) /TP ( 1 ) /Z ( 1 ) ^0. 5«ARE A ( N) /TP(N) /Z CN) 

DO 92 J=2,NJ 

GAS8=GA£R^AREA( J )/TP( J)/Z ( J) 

GASB=GASR>«'C4«P 

GSR ATE = (GAS8-GASA )/0T 
GAS=GASB-GSTART 

GASCHK=GASCHK^0.25^(V (2)<"VH0LD) / (C8 /ARE A ( 2 )) *0T ♦ ( PHCLD/ ZH0L0/T( 2)^ 
1P/Z(2 )/TP( 2) I 


STEP- 14“ 

FIND HEAT FLOW RATE TO WALL AND TOTAL HEAT ADCcD TO WALL 

oo=c9( i)«cpw( n*( Twp( i)“TW( in 

DO 93 J=2,N 

DQ=00+RAO( J I^'SORTI I .O^-ORO X( J ) ♦♦2 ) ♦€ PW (J)*(TWP ( J)-TW (Jn*C9( J) 

0=Q +00 


WRITE OUT RESULTS 
ITtRl^ITEP !♦! 

I F ( ITER 1-IWR ITE) 95,94,94 

CALL WRITE2 

ITERl^C 

IF (TIWE.GE.eNDRHP) CALL WRITE2 
IF (T IME.GE^ENDRMP) RAMP==1000. 


STEP-15- 

CHECK TO SEE IF ENO TIME HAS BEEN REACHED 
IF (T IVE-ENDTIM ) 96,101,101 


ENO TIME NOT REACHED - PREPARE FOR ANOTHER STEP 

IF (T IME .GT.PMl ) DT^O.l 
IF (TIME.GT.PM2) DT=0.2 
IF (TIME.GT.PM3) 0T=0.5 
IF (T IME.GT.PM4) OT*1.0 
DETY=OT 
T IMY=T IME 
TYM E=TIME+OT 
TEIM=TYME + . 1 

IF (TEIM.GT.ENDTIM) TYME^ENOTIM 

CALL INTfcPP ( FLOW,TIME3,MFLOW,TYMd,FLOWP) 

CHL-CK^FLOWP 
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VP*FLOUP/AREA(N ) 

IF (TYMt.EO.ENOTIM) D T=ENOTI M-TI ME 
IF (CFECK.EO.O. I GO TO <57 
DT=OX/(\((N )»VP)*2. 

97 TIME=TIME*OT 

C 

on 98 J=1*N 
T( J ) = TP( J ) 

98 TK( JI=TWP(J) 

GASA=GASB 

Lorw=c 

C 

IF (CHECK-C.) 99f99ilOC 

99 NJ=N-l 
GO TO AO 

100 NJ=N 
N=N n 
NP=N*1 

OROXIN » = (RAO(NP )-RAO( N)> / ( X (NP >-X(N) ) 

DROX(NPI=CROX(N» 

GO TO ^0 
C 
C 

C STEP- 16- 

C END TIME EXCEEDED - INTERPClATc CONDITIONS AT END TIME 

C 

101 R AT IO = (ENDT1M-TIME+OT) /OT 
T IME=ENCTIM 

C 

on 102 J=lfN 

TP( J)=TI J »+RATIO* (TP( J)-T( jn 

102 TWPU l=TW(J»<-PATIO*ITWP(JI-TWUl J 

C 

Q=Q- DO ♦RATIO* DQ 

GAS=GAS-GASB+GASA ♦RATIO*(GASB-GASA) 

GASCHK»CASCHK-GASB+GASA ♦RATIO* IGA SB-GASAI 
C 

C STEP- 17- 

C 

CALL WRITE2 

GO TO 1 
C 

103 WRITE (6,111) 

WRITE (6,112) (V( J),J=1,N,10) 

WRITE (6,113) ISTAR 
GO TO 62 

lOA WRITE (6,11A) ERROR 
WRITE (6,110) 

WRITE (6,112) (TP(J ), J=1,N,10) 

WRITE (6,115) TIME 
WRITE (6,111) 

WRITE (6,112) CV( J) ,J=1,N,10) 

GO TO 89 
C 

C 
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FCPMAT STATEMENTS 


C 
C 
r. 

C 

10‘S FOfiMAT (1HK,7 1H SOLUTION TO CUBIC BY NE WTQN- RH APS ON RECUIRES MORE 
1 THAN 10 ITERATIONS 1 

106 format (1HI,3CX,29H TANK PRfc SSURI Z A TI ON PRQGRAM/IHJ) 

107 format (BOH 

1 ) 

lOa format (IHK,32H ERROR GREATER THAN ACCEPTABLE 1 

109 FORMAT (lFL»5Xfl3H GAS T(R) = F6.1,5H J= 13 ,9H TIME = F5,1|9H 
1SEC0NDS.9H TEST = E13.5t2X,40H COMPUTED TEMP. GREATER THAN BOUND 
2APY 1 

110 FORMAT (1HL*10X,«3H ULLAGE GAS TEMP PROFILE IN PROGRAM UNITS /IX) 

111 format ( IHLt lOXt 108H INSTANTANEOUS ULLAGE GAS FLORATES FROM TOP TO 
1 INTERFACE-NEGATIVE VALUES ARE UNSTABLE-RATES IN PROGRAM UNTS / IX » 

1 12 FORMAT ( 1P8615.7) 

113 FORMAT (lHLtlOX,15H ITERATIONS = 13) 

114 FORMAT (IHL,10X,18H PERCENT ERROR = F5.l,10H PERCENT /IX) 

115 FORMAT (1HL»10X,17H TIME IN RAMP = F5.1,10H SECONDS /IX) 

116 FORMAT (1HL,10X,43H ULLAGE WALL TEMP PROFILE - PROGRAM UNITS /IX) 

117 FORMAT (lHLt2X,17H GAS SUPPLIED = F7.3,5H LBS»2X,14H GAS CHECK 
1= F7.3,5H L8S»2X,14H RAMP TIME = F7.2.9H SECONDS) 

118 FORMAT <1HK,2X,13H GAS FLOW = E13.5.9H LBS/ SECf2X ,15H ITERATION 
IS = I3,2X,17H HEAT TO WALL = E13.5t5H BTU) 

119 FORMAT ( 1FK,50H ONE OR MORE NEGATIVE GAS VELOCITIES ARE COMPUTED ) 
1015 FORMAT (1HK.23H INITIAL ULLAGE GAS = F6.3,5H LBS) 

1050 FORMAT (IFK, 25H INITIAL HEAT TO WALL = F7.1, 5H BTU) 

2015 FORMAT (1FK,23H INITIAL ULLAGE GAS = F6.3,10H KILOGRAM) 

2050 FORMAT (IFK, 25H INITIAL HEAT TO WALL = Ell. 4, 7H JOULE) 

END 
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SMUPATION P VS TI Nt ♦ St CCNOS 
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APPENDIX D 


SUBROUTINES - COMMON TO BOTH PROGRAMS 


SUfiRfltTINE SINTS 
D IRFNS ICN 

1 7GASI30* .TIME1C30» ,PCAT A (30 ) t T IME2 I30> » 

? ELnw(30).TIME3(30) , TSAT ( 30) ♦ T I ME4( 30> t 

3 TBULK(30>*T (ME5( 301 .tiCUTOI 30 > . T I ME6 ( 30 ) , 

4 CIND(30)*TIME7(30>*ZAV( 130)* 

5 T I NU( 30 ), 0 I ST 1 ( 30 ). TG I N0( 30 ). 01 ST2( 30)* 

6 X( 130)*TU50).TP(130>,7W(130).TWP( 150) ,V(i50). 

7 CPI 130>»CPW(150),H(150)*2(150) »21( 150) *Z2(I50)* 

8 RAOl 150) *AREA( 150) .VRAD( 150) «XRAO( 150) *ORDX( 150) 
f 

niPENS ICN 

1 C5( 150).C9( 150). YTICK (50 ) • XXT ( 50 ) * T1 CK (1 50 ) 

. TmHK( 30) *C(ST7(30)« ThCU(30) *DI ST8 ( 30 >. TS TNI 3 0) ,0 1 ST9 1 30 ) » 
TMAL(30)*OIS11(30)*CP8K(150) .TtaBI 150) »TB( 150) »CPCU(i50)» 

4 TWC( 150)*TC(i50).CPST(150).TST(150) *TZ( 150) *CPAL( 150) * 

5 TALSI 150) ,TBCI 150).TTT( 150) ,PRAM( 150) 

r 

CCHMON 

1 X.T,TP*TM.T)«P»V.CP*CPB*T0AS*TlMEl*MTGAS*PCATA.TIME2.)tP0ATA. 

2 f in*<.TIM£3.MfLn)..7SAT,7IHE4,MTSAT .TBULK ,T IME5 .MTBULK, OOUTD* 

3 1 IRE6.ROOUT.C1ND.TIME7.ROIN,N.NP.XN.ULLAGE,RAD1US,DO. 

4 SPtoGT.XMOLEC ,GST ART .HCCNST , T I ME .GAS ,0* G ASCHK , GSRATE , DT . 

5 ZCCNST.h.HMULT.hEXP.YKAO.XRAO.MKAO.TNKteT.CUTR .OIN.UNUMS* 

6 PTxIN.RTGlK.MTICK .MTBAK.MTCU.MTSS.MTAL. 

I TWIND.CIST1.TGIND.DIS72.YTICK.XXT .TW8K.DI ST7. 

8 TWAL.DIS11.TWCU.D1STB.TSTN.0IST9 

C 

COMMON 

1 C ARAD. ADIEU. SPWSS .0 1 A .KCONST . T ADD. AB . AlC. A3S . A5B. WTB. NT 1C. 

2 W13S.KT5B 

C WEEN USING THIS SUBROUTINE EUR THE PRESSURIZATION PROGRAM - USE 

c Tee same dimension and common appearing in tee main program • 

c 

(EIRACIUS .L£. 0.0) GO TO 10 
RADIOS = RAClUS/0.3048 
1C ULLAGE = ULLAGE/0. 304B 
CARaO = CARAD/0.3048 
HCONST = ECCNST/S.lSPHeOS/l.a 
ADIEU = AGIEU/.C9290304 
SPbGT = SPWGT/ 16.018463 
SPbSS = SPwSS/16. 018463 
UIA = CIA/0.3048 
kCCNSl = RCONST/0.304B 
TACC = TADC *1.8 
A8 = AB/.C9290304 
Air = AlC/ .C9290304 
A3S = A3S/.C929C304 
A58 = A56/.C929C304 
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viTf = hTB/. 45359237 
WTIC = hT IC/. 45359237 
wT3S = kT 3S/. 45359237 
t*T5fi s «T5B/. 45359237 
iBtMTGAS .LT. UGC TO 102 
DC 101 J=l.f'TGAS 

101 TGASIJ) = TGASiJl * 1.8 

102 IFIMPLATA .LT. DGO TO 106 
nr 104 J=1.PPDATA 

104 PCATAiJ) = POATAt J)/47.680258 * 0. 1U00000E«^07 
106 IPIPFLOM .LT. DGO TO 109 
00 1C8 J-l.PFLOM 
1C8 FLOM(J) = FLOM(J)/. 3048**3 

109 IFiHTSAT .LT. 1> GO TO 112 
on 110 J»1.PTSAT 

110 TSAT(J) = TSAT(J) ♦ 1.8 
112 IFIMTBULK .LT. DGO TO 116 

00 114 J^l.PTBULK 
114 TBOLK(J) s TBULK(J) * 1.8 

116 IFIMCCLT .LT. DGO TO 120 
00 117 J=1.HU0UT 

117 UOUTO(J) = COUTDl J)/11348.931 

120 IFIMOIN .LT. DGO TO 130 
no 121 J-l.HUIN 

121 OINO(J) - CINOt Jl/1054.3503 * .3048 

130 IFlMThlN .LT. DGO TO 134 
00 131 Jsl.PTWiN 
TmINO(J)- T*INO(Jl*l.a 

131 niSTKJI- OISTKJ 1/0.3048 

134 IFIMTGIN .LT. DGO TO 138 
no 135 Js l.MTGIN 
TGlNO(J)s TGINC(J)*1.8 

135 0IST2(J)= CIST2U 1/0.3048 

138 IFIMRAD *LT. DGO TO 141 
no 139 J^l.PPAO 
YPA0(J1= yPAO( J 1/0.3048 

139 XRAOIJD XRAOl Jl/0.3048 

141 IFIPTICK .LT. DGO TO 144 
no 142 J=1.MTICK 
yTlCKUl= yTlCKl Jl/0.3048 

142 XXTUi- XXT( J1/C.3048 

144 IFIPTBAK .LT. DGO TO 147 
00 145 J-l.HTBAK 
TtNCKlJls TkBK(Jl*1.8 

145 0IST7UD 0IST71 Jl/0.3048 

147 IFIPTCU .LT. DGO TO 150 
OC 148 J-l.MTCU 
rwruijD TMru(ji*i.8 

148 0IST8UD 01ST81J 1/0.3048 
15C IFIPTSS .LT. DGO TO 153 

on 151 J=1.PTSS 
TSTMJ1= TSTN(J1*1.8 
151 0IST9U> = DIST9( Jl/0.3048 

153 IFIMTAL .LT. DGO TO 156 
no 154 J=1.PTAL 
TRAKjls T8AL(J1*1.8 

154 OISIIUDOISIK Jl/0.3048 
156 CONTINUE 


RETURN 



n r> r r- 


C fCH THE PRESSOR I/ATICN PROGRAM - REMOVE ALL STATEMENTS FROM 144 

r THROUGH 154 - CHANGE STATEMENT NUMBER 156 TO 144 • 

C 

ENT 


SUEftfUTCNE SPHEATITt TW.CPfCPHtXMOLEO 
C 

D I MENS ICN CPTERPi 15I,CPMALL( 15) *CTEMP( 15) tCGAS( 15 ) 

C 

DATA FOR SP HT OF ALUM ALLOV 2^1S-T852 

< THFRMOPHYSICAL PROPERTIES RESEARCH CEN TER-PUROUE UNIV#) 

DATA (CPTEMP( J ) tCPWALL (J ) • 1 1 15 ) / 36* # .0020 1 45 • • .0042 1 54, » .0074# 72 

1. • .016 U9C...03 30 « 108. # .0513# 126. « .0692# 144. #. 0837 #162 •«. 0994# 180. 

2. . !b2C#2 70. .. 16CC.315. .. I 72C#360. # .1830 #432# #• 202 9 # 540 . # . 2080/ 


C 

c 


c 


f 


c 


f 


DATA FOR SPECIFIC HEAT QF HYDROGEN AT 50 PS lA (N8S) 

DATA (CTEMP(J)#CGAS( J) #J=si # 15) /45. 5# 3.65 #48. 8# 3. 40 #50. 0 , 3 . 30 • 55 . 0 
1 #3.05 #56. 0.2.92# 60. 0#2 .65# 65.0# 2 .76 #70.0 #2 .6 5# 80.0 #2.o2# 100. 0# 
22.6C#200.0.2.72#25U.O#2.d8#300«0#3.0b#350.0#3.28«450.0#3.42/ 


IF ITfci-CPTEMP( U) 1#1#3 
I CPW = CPWALLI 1 ) 

GO TC 9 

3 DC 7 1^2.15 

IF (Tk-CPTFMPU>) 5#5#7 

5 CPk = CPbALUI-l) 3 4TW •CPTEMPn-l))/(CPTEMPU)-CPTEMPl I-l) ) 
1 ♦(CP WALL < I )-CPWALL ( I-l ) ) 

GC TO 9 


7 CONTINUE 

CPk = CPWALU 15) 

9 IF (XMOLEC-3.0) 13#13#ll 
11 CP = 1.24 
RETURN 


13 IE (T-CTEMP(D) 15#15.17 
J 5 CP =: CGAS( ) ) 

RETURN 

17 or 21 J-2# 15 

IF (T-CTEMP( Jl > 19# 19,21 

19 CP = CGAS(J-U^(T-CTEMP( J-1) )/(CTEMP( J)-CT£MP( J-1) )♦ 
1 (CGAS(J)-CGAS(J-l)) 

RETURN 


r 

2 I CCNT IMJE 
C 

CP * CGAS( 15) 
RETURN 


r 

f 

f 


ENC 
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r 

£ 


SUBROUTINE MRITE 1 


ni»*ENS ION 

1 TGAS(30l«TIMEU30)*PCATA(30).TIME2Ii0) • 

ELOt«I30).TIHE3(30)*TSAT(30) »TIME4(33I * 

TBULK(30).T IME5 ( 30 > .OGUTOI 30 } .T IME6 ( 30 )» 
CIN0I30)«TIME7(30)«ZAW( 1501* 
TMlNC(30).0ISTlt30)tTGIN0(30).0IST2(30)» 

XI 150).T( iSO)*TP(150).TMU30>.TMP(150) »V(150)> 

CPI 150)«CPta(lS0)«Hll50>»Z(150)*Zi(130) *Z2 1 150) » 

RAOl 150) « AREAI 150) •YRAOl 150) •XRAD(150) tOROXI 150) 

UIPENSICN 

C5I 150)*C9( 150)«VTICKI50)*XXT(50) .TICK (150) 
.TmBKI30)*OIST7(30),TaCU( 30)«OIST8(30) .TSTN(30)*0IST9(30) t 
TtaAL(30).0(Sll(30).CPBK( 150)*TWB( 150) »TB(150) .CPCUI150)* 
TtoCI 150),TC( 150),CPST( i50),TST(l50) «TZ( 150) .CPALl 150) , 

TALSI 150) ,T8C( 150 ) . TTT ( 1 50) .PRAHI 150) 

CCPMCX 

X«T.TP*TM«TWP«VtCP«CPK.TGAS*TIHEl.MTGAStP0ATA.TIME2.MPDATA. 
FL0ta.T(HE3.HFL0k*TSAT*TIME4*MTSAT.TBULK,TIME5»MTBULK*00UT0, 
1 iME6*MOOUT*OINO«TIME7.PUlN*N*NP.XN*ULLAGE,RAOIUS.UU* 
SPMGT.XMDLEC.GSTART.HCONST.TIME.GAS.OfGASCHKtGSRATEfOT* 

5 ZCCNST.H.HPULT.HEXPt YRAO.XRAO.MRAOt TNKMT.CCTR.OIN.UNUMS. 

6 PTmIK.PTGIN.MTICK.MTBAK.HTCU.MTSS.MTAL* 

7 TWIN0.DISTl.TGIN0.0IST2«yTICK,XXT.T)»8K.CIST7, 

e TtaAI. .nisi 1.TMCU.0IST8.TSTN.0IST9 

C 

CCHMON 

1 CARAO.AOIFU* SPMSS *01 A .RCONST * T ADO . A £ » A1 C. A3S . A5B. IkTB. MTIC . 

2 t»l3S.KT5B 
f 

C WHEN USING THIS SUBROUTINE FCR THE PRfcSSURIZAT ION PROGRAM - USE 

C. TFE SAME DIMENSION AND CCMMON APPEARING IN THE MAIN PROGRAM . 

f 
f 

MRITF (6.1007) 
f 

1007 FORMAT I IHL . 12H INPUT DATA) 

C 

N = XN 

UX = ULLAGE/I XN-1.0) 
nC 30 JsP.N 
XJTICK = J-l 

30 X(.J> - CX « XJTICK 

X( 1 ) :: 0. 

no 31 J*1.N 

CALL INTFRP ( TM IN C.U ISTl . MTW IN .X( J ) . Tm( J ) ) 

Jl CAi L INTERP (TGIN0.01ST2.MTGIN .X( J).T( J) ) 
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IFUJ^Ll^*S i 100# lOl t 100 

iUl rtRire (btlOiOl XNfbLLAGE^RAOlUSf SPiiGT,XMOL£C 

r 

1010 FCRHAT nHl,9X.3^h INITIAL NUf^BER OF NETPOINTS = F3.0/ 

1 10X,ISIH INITIAL ULLAGE = F6*3#t>H FEET/ 

/ 10X.16H TANK KAOIUS = F6.3.6H FEET/ 

3 lOXtBOH TANK UALL SPECIFIC WEIGHT = F6.l,16H LBS/CUtilC FOOT/ 

A lOXtPlH MOLECULAR WEIGHT = F3.1I 
WRITE I6^10C8) 

C 

iOOH FORMAT IlHlf9X.4CH GAS T6MP#DEG R VS T I ME • SEC ONDS / IX 1 

r 

WRITE (6«1009I (TGAS(JI« TIMElIJlt J - IfMTGAS) 

f 

1C09 FCRMAT IF20.1* E23.ll 
C 

WRITE (6.10301 

r 

1030 FORMAT (lhL.9X.47H PRESSURE.LeS/SO FT VS TIME. SECONDS 

1 /IXI 

r 

WRITE (6. 10091 CPCATAI J1 .TIM E2 1 Jl.J^l.MPOATAl 

C 

WRITE (6.10311 

r 

IU31 FCPMATI IHL .9X.44F FLOW RATE. CU FT/SEC VS TIME. SECONDS / 

11X1 

c 

WRITE (6.10321 ( FLOW! Jl .TIME3IJ1 .J=1.MFL0W1 

r 

1032 FORMAT ( F20 .4.F25 . 1 1 

f 

WRITE 16.10331 

1033 FfRMAT (IHL.9X.49H SATURATION TEMP. DEG R VS TIME. SECONDS 

1 / 1 X 1 

C 

WRITE (6. 10091 (TSATI J1«TIME4( Jl .J-l.MTSATl 
f 

WRITE (t. 10341 
C 

1034 FORMAT ilFL.9X.44H BULK TEMP. DEG R VS TIME. SECONDS /I 

iX 1 

r 

WRITE I6.1C091 ( TRULKU 1 .TIMESIJl .J-l.MTBULKl 

f 

WRITE (6.10351 
f 

1035 FORMAT (1HI.9X.63H CUTSIOE HEATING RATE.8TU/SC FT-SEC VS 

ITIME. SECONDS /IX) 

C 

WRITE 16.10321 (00UTD4 Jl .TIME6( J) .J=l.MuOUTl 

C 

WRITE (6.1036) 

r 

1036 FORMAT (IHL.9X.62H INSIDE HEATING RA T£ . BTU/ LI NEAR FT-SEC VS 
ITIME^SECONCS / IX 1 
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f 

UKITE (6,103?) (CINOt J),T(ME7(J) ,J=1,MU1N) 

C 

mRITE (6.10)3) 

C 

1013 FCEMAT UHL,9X*3SH INITIAL TEMPERATURE-OEG. RANKINE /IHK, 

1 9X.14H X CUUR0INATE,5X. 16H WALL TEMPERATURE ,5X * 

2 17H GAS TEMPERATURE/IX) 

f. 

WRITE (6.101A) (X( J).TW( J) ,T( J), Jsl.N) 

C 

1014 FORMAT (E1G.2, F22.1, F23.1) 

f 

WRITE (6.20S0) 

20^0 FORMAT! IHL,9X,48F TANK RADIUS (FEET) VS AXIAL DISTANCE (FEET) ) 
WRITE (6,2051) ( YRAOI J ) , XRAD ( J ) . J^l.MRAO) 

2051 F0RMA1 ( F20 .3.F25 .2 ) 

r 

IF (FCCNST) 5.5,7 

f. 

5 WRITE (6.1039) 

C. 

1039 FORMAT (IHK, 46H HEAT TRANSFER COEFFICIENT WILL BE COMPUTED ) 

WRITE (6,1051) HMULT, HEXP 

1051 FORMAT (IHJ.lOH HMULT « F6.4,9H HEXP = F6.4) 

GO TO 9 
C 

7 HCCNST = 3600.4HCCNST 
WRITE (6,1040) HCCNST 
C 

1040 FORMAT (IHK, 33H CONSTANT HEAT TRANSFER COEFF = F6.2.23H BTU/SO 
IFT-HOUfi-CFG R ) 

C 

HCCNST = HCONST/3600. 

60 TO 9 

100 WRITE (6,4010) XN, ULLAGE, RADIUS, SPWGT.XMOLEC 

C 

4010 FORMAT (1HL.9X.32H INITIAL NUMBER OF NETPOINTS - F3.0/ 

1 10X.19H INITIAL ULLAGE * F7.3.7H METER/ 

2 10X.16H TANK RADIUS = F6.3.7H METER/ 

3 10X,30H TANK WALL SPECIFIC WEIGHT = F6.1.16H KG/CUBIC METER/ 

4 10X.21H MOLECULAR WEIGHT = F3.1) 

WRITE (6,4008) 

0 

4008 FORMAT (1HL,9X,4CH GAS TEMP.OEG K WS T I ME , SECONDS / IX ) 

WRITE (6,4009) (TGAS(J), TIMEKJ). J = l.MTGAS) 

4009 FORMAT (F20.1, F25.1) 

0 

WRITE (6,4030) 

4030 FORMAT (1HI,9X,44H PRESSURE , M-NEWTON/SU M VS TI ME , SECONDS / IX ) 
WRITE (6,4009) ( POATA( J ) ,T IME2 ( U > , I , MPDA TA > 

WRITE (6.4031) 

^.031 FORMAT ( IHL ,9X,44H FLOW RATE, CU M/SEC VS TIME, SECONDS / 

1 IX ) 

WRITE (6.4032) ( FLOW( J) , T I ME 3 ( J ) , J=1 , MFLOw) 

'^032 FORMAT ( F20 .4, F25 . 1 ) 

wRITF (6.4033) 
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^.0^^ FflJ^MAT UHL.9X.49H SATUKATIGN TEMP»0EG K VS TIME, SECONDS 

I / IX I 

mPITE (6.4CC9) (TSATI J),TIME4(J),J=1,MTSATI 
MhITE (6,4034) 

4034 ECRMAT (IHL,9X.44H BULK TEPP,DeG K VS TIME, SECONDS /I 

IX ) 

WRITE I6,400S) (TBULK(J),TIME5( J> ,J=l*MTdULK) 

WRITE (6,4035) 

4035 FORMAT <1HL,9X,63H OUTSIDE HEATING RATE , JOULE /SO M-SEC VS 

IT IRE, SECONDS /IX) 

WRITE (6,4032) (OOUTD(J) ,TIME6( J) ,J=1,M00UT) 

WRITE (6,4036) 

403b FORMAT (1M,9X,62H INSIDE HEATING RATE, JOULE /METER-SEC VS 

IT IRF,SECONCS /IX) 

mRITE (6,4032) (UIND( J),T1ME7(J) ,J=1,M0IN) 

WRITE (6.4013) 

4013 FORMAT (IHL,9X,34H INITIAL TEMPE RATURE-DEG. KELVIN /IHK, 

1 9X.14H X COORDINATE, 5X, 18H WALL TEMPERATURE ,5X, 

2 17H (,AS TEMPERATURE/ IX) 

C 

WRITE (6.4014) ( X( J),TU( J) ,T( J), J-l.N) 

4014 FORMAT (F19.2, F22.1. F23.1) 

IF (HCONST) 50,5C,70 

60 WRITE (6.4039) 

4039 FORMAT I IHK. 46H HEAT TRANSFER COEFFICIENT WILL BE COMPUTED ) 

WRITE (6.4051) HMULT, HEXP 

( 

4061 FORMAT (IHJ.lOH HMULT = F6.4.9H HEXP = F6.4) 

GO TO 90 

70 WRITE (6,4040) 

4040 FORMAT (IHK, 33H CONSTANT HEAT TRANSFER COEFF = F6.2.24H WATTS/S 
lU-METER-DEG K ) 

SO f.CNTINUE 

WRITE (6. 4050) 

-4060 FORMAT (1HI,9X,48H TANK RADIUS (METER) VS AXIAL DISTANCE (METER) 
I ) 

WRITE (6*2051) (VRAD( J).XRAC( J), J-1,MRA0) 

C 

9 CfKTINUE 
C 

IF (/CONST) 13,13,15 
13 wRITF (6.1042) 
f 

1042 FORMAT (IHK, 25H A REAL GAS IS ASSUMED ) 

GO TO 17 

r 

15 WRITE (6.1041) 
f 

1041 FORMAT (IHK, 27H AN IDEAL GAS IS ASSUMED ) 

f 

17 CONTINUE 
f 

IF (XMCLEC - 3.0) 22,22,23 
22 WRITE (6,1011) 
f 
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lOl L FORMAT (IhK, 28F 
C 

(;c Tf 2 5 
r 

2^ wRITF (6,1012) 
f 

1012 FfiRMAT ilHKt 26F 
f 

25 CCMUUF 
RETURN 
C 

FNC 


TFE PRESSURANT IS HYOROGEN/lHl *9 X»9H RESULTS) 


THE PRESSURANT IS HEL I UM/lHl ,9X t 9H RESULTS) 
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-> n n 


SOBRCUTINE COMPRS (TtP*Z2t2U/2^XMCtEC) 
C 

c 

Oi^ENSCCN PX(17)tTX(20)fZ(20fl7> 

C 

I IP ( 3 . 0 -XMCLEC) 3 f 5.5 

H // = 1 .0 

21 = 1*0 
72 - 1*0 
C 

RtTURR 

C 

5 CCMINUE 


f 


COMPRESSI eiLlTY 
NBS TN 120A) 


DATA IZITfP) VALUES COMPUTED FROM 


DATA 


OF 


OATAIPXIJ 17)/ 


0ATA(TX<J)*J=^1»20>/ 


8640^ • 10080* 9 
20160*«23040*» 
34*« 38*« 42* 
66*9 70*9 /4* 

I20*9l4099lb0* 
0ATAUZIJ9K >9J^l9 20>9K=^l9 9)/ 

1 *9153 9 *9345* *9482* *9 577 

2 * 9825 ** 98389 * 98499*9858 

1 •8766* *9074**92769*9411 

2 *9782**9 801**9818* *9832 

1 * 87189 * 87189 * 90129*9212 

2 *97279*9754**9777**9796 

1 *8036* *8036**848 7* *8819 

2 *962 9* *9670**9 704* *9 733 

1 *80 92**8092* *80929*8391 

2 *9529**9585**9631**9670 

1 *729 7**72979*7297** 7916 

2 *9429 9*9499* *9558* *96C6 

1 *739 2**7 392* *7392* * 7392 

2 *9327**94139*9484**9543 

1 *682 7**6827* *6 827* *682 7 

2 *9224**9326**9410* *9479 

1 *6281**628l**6281**6281 

2 *9120**9239**933 5* *9415 
CATA( (ZU«K >9 J=)*20)*K = 10* 17) 

1 *6 593**6 593* *6593* *6 59 3 

2 *9014**9150* *9260** 9351 

1 *6735**6735**6735**6735 

2 .8906 9*906 1 9*91 8 5* *9287 

1 *609 8**6098**6C9 8* *6C96 

2 *8699**8876**9031* *9157 

1 *4689**4689**4689**4689 


2880*9 4320*9 5760 ** 7200 ** 
11520**1296 0 * * 14400 ** 17280 * * 
25920 * *28 600 * * 36000 */ 

♦ 46 ** 50 ** 54 ** 58 ** 62*9 

• 78 ** 82*9 86*9 90 ** 100*9 
* 200 */ 


*9644 * 
* 9866 * 
* 9513 * 
. 9843 * 
. 9355 * 
* 9312 * 
* 9049 * 
* 9757 * 
* 8727 * 
* 9702 * 
* 8380 * 

* 9<^4 7 * 

• 6002 * 
* 9592 * 
• 7 586 * 
* 9537 * 
* 7114 * 
* 9482 * 

* 6593 * 
* 9427 * 
• 6735 ♦ 
* 9372 * 
* 6098 * 
* 9262 * 
*4689 * 


9695 * 
9882 * 
9589 * 
9866 * 
9461 , 
9643 * 
9217 * 
9803 * 
8963 * 
9164 * 
8697 * 
9725 * 
8416 * 
9686 * 
8118 * 
9647 * 
7 79 8 , 
9608 * 

1450 * 

9569 * 

7063 * 

9531 * 

6098 * 

9452 , 

4689 * 


9734 * 
9900 , 
964 7 * 
. 9893 * 
. 9542 * 
. 9880 * 
. 9343 # 
. 9659 * 
. 9139 * 
. 9839 * 
. 8928 * 
9819 * 
. 8708 * 
. 9798 * 
8479 * 
. 9778 , 
. 8240 * 
. 9758 * 


* 9764 * 

* 9910 * 

* 9693 , 

* 9908 * 

* 9605 , 

* 9899 * 

* 9441 , 

* 9689 * 

* 9273 * 

* 9678 * 

* 9102 * 

* 9868 * 

* 8926 * 

* 9858 * 

* 8744 * 

* 9848 , 

* 8557 * 

* 9837 * 


. 9789 , 

* 9916 * 

* 9729 , 

* 9916 , 

* 9654 * 

* 9911 * 

. 9518 , 

* 9906 * 

* 9378 * 

* 9901 , 

* 9236 , 

* 9896 * 

* 9092 * 

* 9892 * 

• 8945 * 

* 9887 * 

* 8794 * 

* 9883 * 


* 9809 * 

. 9922 , 

* 9758 * 

. 9924 , 

* 9695 * 

* 9922 * 

. 9579 , 

. 9923 , 

* 9462 , 

. 9923 , 

* 9343 * 

. 9924 * 

. 9222 * 

• 9925 * 

* 9100 * 

. 9926 * 

* 8975 * 

* 9926 / 


* 7989 ** 8364 * 

* 9739 ** 9828 * 

* 7723 ** 8163 * 

. 9719 ** 9818 * 

*7136**7733* 

* 9680 ** 9798 , 

. 6429 ** 7256 * 


*8641**8848f 
.9878*. 9927* 
*8484**8718* 
*9874**9928, 
• 8145* *8462* 
*9866**9932* 
*7775**8204* 


2 

1 

2 

1 

2 

1 

2 
1 
2 


• 8492 ** 8671 ** 8866 * * 9 C 24 
* 3298 *. 3298 . * 3298* *3298 
* 8233 *. 8530 *. 8 734 * *86 74 

• B 632 ** 3632 ** 3632 * *3632 

• 7952 * * 8401 ** 8603 * *86 55 
. 5120*. 5 1 20 . . 5120 * *5120 
. 7792 ** 8123 * . 6412**8654 

• 6427 **6427 ** 6427 * *6 42 7 
. 7150 ** 7661 ,* 8028 * *8314 


* 9 l 52 « • 9372 * * 96^ 1 * *9 Z 80 * 
* 3298 ** 3298 * * 5332 * *677 4 * 
* 9042 ** 9288 * *9 602 **9 761 * 
* 3632 * * 3632 * * 3632 * * 6260 * 
* 89 32 ** 9188 , * 9563 * *9 744 * 
* 5120**512 0 * * 5120 ** 5120 * 
. 8823 * *902 7* .9 524 * *9 726 * 
*642 7 ** 6427* .642 7 * *6 427 * 
* 8553 * * 9001 * * 9414 * . 9684 * 


9858 ** 9937 , 
7476 ** 7882 * 
9850 ** 9942 * 
7232 ** 7520 * 
9843 ** 9947 * 
6 ^ 94 ** 7327 * 
9836 * * 9953 , 
6427 ** 6427 * 
9821 ** 9960 / 
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IF (TX(20>-T) 




11 = UO 
1 \ = 1#0 

12 = 1.0 

RFTUftN 

or 17 j=i.i^ 

JJ = ao - J 

IF <TX«JJ>-T> 15f 15fl7 

NT = JJ 
GO TO 19 

CCKT IMJE 
GO TO 11 

IF<P-PX( 1 > > 11.23*23 

00 27 J=2.17 

IF (P-PXIJl) 25.25.27 

NP J-1 
GO TO 31 

CCM INUE 
GO TO 11 

Z A = Z INT.NPIX P-PX<NP)I / C FX ( NF^-1 >-PX ( NP ) I ♦< Z( NT .NP^l ) -Z (NT . NP ) ) 

ZR = Z<NT41.NP) + IP-PX(NP))/CPX(NP+1)-PX(NP)I4^ 

1 I ZINT+l.NP^l )-Z(NT+l.NPll 

7Z = ZA 4 IT-TX(NT))/(TX(NT+l)-TX(NTl)*(Za-ZA) 

U/CT ^ (ZB-2A)/<TX(NT^I)-TX(NT)> 

Z1 = ZZ^T^CZDT 

OZCP ^ (Z<NT.NP+l)-Z(NT.NP)+(T-TX(NTI)/(TXiNT+ll-TX(NTl)* 

1 ( Z iNT^i.NP^l )-ZINT+l.NPl-ZINT.NP+l)^Z(NT,NPI) ) / ( PX ( NP+1 ) -PXI NP )) 

17 = ZZ-P*CZOP 


RFTiJPX 



SUJ^Rf UTlNF hCOEFF ( T f TTT t P • XL t Ht PR AK • ZCCNST • LL T • HEXP tXMOLEC) 
r 
( 

Or^ENSCGN VIY(2,7),VIX(2*7>*C0Y42»7) •COX(2t7l 

r 

C 

C FYORGGEN viscosity data ( (LBM SEC) /( so FT)) 

C ICATA FROM HILSENRATH* ET AL ) 

OATA(VlY{ l9J)«VlX(ltJ)fJ^l«7) / •l064E“7fi6» t*6Q0S E** 7 # 1 08« ? 

1 1.001E-7*216. • 1.324E-7*324«, 1.65>4E-7,45 0# »1 .9 54E-7,576. • 

2 /.195E-7*<;E4./ 
f 

f 

C hydrogen thermal conductivity DATA ( BTU/ ( SEC > ( FT ) ( OEG R/FT)) 

C ICATA FROM HlLSENRATH* ET AL ) 

CATAICOYI l^J)^CCXI !• J) tJ=i,7) / .il87E-5tl8. ».6780E-5» 108. » 

1 1.26 8 E- 5 . 216 .. l.87ie-5f 324. , 2. 505E-5 .41>0. » 3. 07 5E-5t 3 76. t 

2 3.520E-5.Z84./ 

C 

c 

C FELIIJM VISCOSITY DATA (LB-SEC/SO FT) (VANCE AND DUKE) 

r 

OAT A (VIY(2fJ).VlX(2.«i).J = lf7)/ .376E“7tl0.f .668E— 7 1 30. t i. 002E*7 1 
I 60.. 1 .420E-7f 100..2.257E-7.2C0. •3.528E-7.400..6.075E-7.800./ 
f 
C 

C HELIUM THERMAL CONO DATA I B TL/ SEC-FT-OEG R) (VANCE AND DUKE) 

C 

CAT A (C0Y(2.J ).CCX(2. J).J=1.7)y . 1806E-5 . 10 . t .373E-5 .30. . . 397E-5. 

I 60.. .834E-5. 1CC..1.277E-3,20C. .2.E-3f400..3.445E-5.800./ 

I M = .S^XMGLEC 
R = 1545.4 

C = XM0LEC4P/R 
C 

TAV = 0.5*(T^TTT) 

RHCSO = (C/TAV)#«2 
C 

DC / J=2.7 

IF (TAV-VIXIM. J)) 3.3.7 

C 

3 Vise = VlY(M.J-l)+(TAV-VIX(M.a-l))/( VIX(M.J)-VIX(M.J-i)>* 

1 (VIY(M.J)-VIY(M.J-I)) 

GO TO 9 
C 

7 CCNTINUE 

Vise = VIYIM.7) 

C 

9 on 15 J=2.7 

IF (TAV-COX(M. J)) 11.11.15 

r 

II CCNO =5 COY(M, J-1)+(TAV-C0X(H. J-1))/(C0X(M. J)-CCX(M.J-l))* 

1 (CGY(M. J)-COY(M. J-D) 

GO TO 17 
C 

15 CCNTINUE 

COND C0Y(M.7) 

r 

17 CALL SPHEAT ( T AV . TAV .CP . CP . XMGLEC ) 

r 
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J 


TUIFF= AhS(TTT-T) 

IF (TOIFF-1.0) IS, 23, 23 
IS TOIFF = l.C 
23 IF (ZCCNST) 27,27,25 
25 HF73 * \,onik\l 
C,r TO 2S 

27 call CCMPRSi TAV,P,Z,21,Z2,XI>CLEC) 

BETA = Zl/Z/TAV 

2S XHCLO * Al.OG(XL**3*RFOSO*CP*TOIFF/VlSC/CCNO*BETAI 
XhCLC = HEXP*XHOLD 
XFCLO = FXP(XHOLD) 

H = HPLLT*<CONO/XL)*XHCLO 
PVC - CCNO ♦ ( 1.0/(32.17*VISC))**0,8 
POL = (CP ♦ VISC/CCNO * 32.17)**.333 
PRAM a POL ♦ PVC 

METURA 

ENT 
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St.BKCUTlNE INTEPP ( Y * X t M « AKG* ANS > 


OIHENSICN Y( 100),X(100) 

IF (t^-n 3t5,7 

KETUKN 

ANS = Yd) 

RETURN 

IF (ARG-XI in 5«5*9 

on )i j*2.)» 

IFIARG-X(J)) mil. 13 

AIMS = Y(J-1)*(ARG-X(J-1))/(XIJ)-X(J-1))*(V(J)-Y(J-1>) 
RETURN 

CENT INUE 
ANS > VIM) 


RETURN 



OUTPUT FOR EXPULSION PROGRAM 


StBROUTlNE WRITE 2 


C 

ni»*ENSiON 

1 TGAS(30)«TIME1I30I.PCATA(30}.TIHE2(30I . 

7 FLQta(30)«TIME3(30) tT SAT I 30) »TIME4(30). 

3 TBULK(30)«TfHE5(30).OOUTO(30).TlHE6(30)» 

4 C 1N0430).TIRE7430)*2AV( 150). 

5 TwCNC430)«0(ST1430).T61N0430)*OIST2430)t 

6 X4150)«T4150)*TP4150)*TW4150>.TMP4150) •V4I50)* 

7 CP4 150).CPh4150)«H4150)*24l50)«ZI4150) *224150)* 

8 RA04 150).AREA4150).YRA04 150) «XRAO 4 ISO) *0R0X4150) 

OtPENS ICM 

1 C54150).C94 150)*YTICK4S0)*XXT450)*TICK4150) 

2 *TWBK4 3O)*0IST7 4 30)«T)iCU4 30) *OIST84 30) * TSTN430) *0IST9 430) * 

3 TWAL430) *0IS11430)*CPBK4 150)*TWB4 150) *TB4 150) *CPCU4150)* 

4 TWC4 150)*TC4 150)*CPST4 150)*TST4 150) *TZ4 150) *CPAL4 150) * 

5 TAI.S4 150) *TBC4 150)*TTT4 150) *PRAN4 150) 

C 

CCMMON 

1 X.T*TP*TW*TWP.V*CP,CPW»TaAS«TIMEl*HTGAS*PCATA*TIME2*NP0ATA. 

2 FLaw«TlHE3*MFL0)>*TSAT.TINE4.HTSAT.TBULK*TlME5*MTBULK*U0UT0* 

3 T I PE6.ROaLT.01NO*TlME7«POIN*N*NP.XN .ULLAGE* RADIUS *00* 

4 SPWGT.XMOLEC.GSTART.HCONST*TIME.GAS*0*GASCHK«GSRATE.OT* 

5 ZCONST«H*HPULT«HEXP*YRAO*XRAO.MRAO*TNKMT*COTR*OIN*UNUMS* 

6 PTMlN*PTGIA*MTICK*MTeAK.NTCU*MTSS*MTAL* 

7 lWlNO.OiSTl.TGIK0.01ST2«YTICK«XXT.Tt«BK*OtST7* 

8 TWAL.0IS11*TWCU*0IST8*TSTN.QIST9 
f 

COHMCN 

1 CARA0.A0IFU*SPWSS*01A*RC0NST*TA00*AB*A1C* A3S*A5B*UTB*MT1C* 

2 aT3S*)>T5B 
C 

IF4UKLHS .GT. 0.) GO TO 200 
WRITE 46*1016) T1ME«GAS*0*V42) 
f 

1016 FORMAT 41HL/1HL*9H TIME - F6.1.9H SECCNOS* 5X *17H GAS SUPPLIED 

1 F7.3.5H LHS*5X*17H HEAT TO MALL = FB.1*5H BTU*5X* 

2 14F INLET VEL * F6.4.10H FEET/SEC) 

C 

WRITE 46.1017) COTR.GASCHK.GSRATE.OIN 
C 

1017 FORMAT 41HK.16H HEAT TO LIO < F7.1.5H BTUS*3X*17H CHECK ON GAS 

1 F7.3.5H LBS.3X.13H GAS FLOW = F7.4.9H LBS/SEC* 

2 3X.16H HEAT TC HARD * F6.4.13H BTU/4FT-SEC )) 

C 

WRITE 46.1018) 
f 

1018 FORMAT 41HK.10H X. FEET.2X.11H WALL T4R)*11H GAS T4R)*7X* 

1 lOH X. FEET.2X.11H WALL T4R).11H GAS T4R)*7X*10H X* FEET 

2 2X.11H WALL T4R)*11H GAS T4R)/IX) 

C 

WRITE 46.1020) 4 X4 J ) • TWP 4 J ) . TP4 J ) ♦ J= 1 .NP ) 

1020 FORMAT 4F10.2.F12.1.F12.1.F17.2.F12. 1 .F 12. 1 . FI 7. 2 . F 12. 1.F12.1) 

C 

RETURN 

20C CfNT IMJE 



on /Oi j=i*NP 
X(J)= X(J)*0.3C48 
TP(J)= TP(J)/1.8 
TkiP(J)= TteP(J)/1.8 
201 CCMIMjE 

GAS = GAS * 0.4535S237 

RT = C ♦ 1054.3503 

V(2I- V(2)«0.3048 

COLA = CCTR ♦ 1054.3503 

OeiCK = GASCHK ♦ 0.45359237 

GSRATE = GSRATE ♦ 0.45359237 

UIN= C IN 4 3459.1475 

WRITE 16.300) TIME«GAS*RT*W(2) 

3U0 FORRAT (1HL/1HL.9H TIME » F6.1.9H SECCNOS. 5X .16H GAS SUPPLIED = 
1F7.4.9H KILOGRAM. 1X.17H HEAT TO WALL = E1L.4.6H JOULE. IX. 14H INL 
PET VEL * F6.4.1CH METER/SEC) 

WRITE (6.305) COLA. OBLCK. GSRATE. OIN 

305 FORMAT (1HK.16H FEAT TO LIO - E12.5.2H J.1X.17H CHECK ON GAS - F 
17.3.4H KG.4X.13H GAS FLOW » F7.4.8H KG/SEC. IX. 16H HEAT TO HARO 

2- Ell.A.lOH J/(M-SEO) 

WRITE (6.310) 

31C FORMAT (IHK.lOH X. METER. IX. 13H WALL TEMP-K.14H GAS TEMP-K.3X 
i.lOH X. METER. IX. 13H WALL TEMP>K.14H GAS TEMP-K .3X. lOH X. ME 
2TER.1X.13H WALL TEMP-K.14H GAS TEMP-K/IX) 

WR(TE(6.311) (X(J).TWP(J).TP( J).J=1.NP) 

311 FORMAT (F10.4.F12.1*F12.1.F17.4.Fi2. I . F I . I , F 1 ? . . F12. 1 . FI 2.1 ) 

210 COAT INUE 

00 21 1 J=1.NP 
X(J) 3 X(J)/.3048 
TP(J) ® TP(J) ♦ 1.8 
TWP(J) 3 TWP(J) * 1.8 

211 CCKTIKUE 

GAS • GAS / 0.45359237 
V(2) - V(2)/.3048 
OIN = CIN/3459. 1475 
C 

RETURN 

r 

ENC 
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OUTPUT FOR RAMP PROGRAM 


SUBR0L1 INE WR ITE2 
C 

DIMENSION TGASIlOOJt TIMEKIOOJ, POATA(IOO), TIMEZaOO), FLCW(25), 
I TIME3(25I, TSAT(50I» TIME4(50), TBULK(25), riM65(25), CO!JTD(25lt 
2TIME6(25), QINO(25), TIME7(25), ZAV(250)t TWIND(100>f DISTKIOJ), 
3TGIN0(100), DIST2(lOOI, X(25C), 1(2501 , TP (250 > , TW (250 > , TWP(250) 
4, V(25C), VP( 250), CP(250), CPW(250), H(250), Z(250 ), Zl(250), Z 2( 
5250), RAD(250), AR6A(250), YRAD(250), XRA0(250), DRCX(250), TtST( 2 
650) 

DIMENSION C5( 250), C9(250), YTICK(250) , XXT(250), TICK(250), XR(3) 
I, A0(3) 
r 

common X,T,TP,TW,TWP, V,CP ,CPW,TGAS,TIME1 ,MTGAS,PDATA,TIME2, mpoata, 
IFLOW, TIME3,MFL0W, TSAT , TI ME 4 , MTSA T , TBULK , TI ME5 , MTBU LK , QGUTD , T I ME6, M 
2Q0UT, 0IND,TIME7,MQIN,N,NP ,XN,ULLAGE,RADIUS,TLIO,SPWGT ,XMOLEC, GSTAR 
3T,HCONST,TIME,GAS,Q,GASCHK,GSRATE,DT,ZCONST,H,HMULT,HEXP,YRAD,XRAD 
4, MR AD, CQ,QIN,UNUMS,TWIND,DIST1,MTWI N,TGIND,DI ST2,MTGI N,0IA, 

5TAOO, YT 1CK,XXT,MTICK 
C 

IF(UNUM£ .GT. 0.) GO TO 200 
WRITE (6,1) TIME, GAS, Q, V( 1) 

C 

C 

WRITE (6,2) GASCHK,GSRATE ,V(NP) 

C 

C 

WRITE (6,3) 

C 

C 

WRITE (6,4) (X(J),TWP( J),TP(J) ,J=1,N) 

C 

R ETURN 

200 CONTINUE 

DO 201 J=1,N 
X(J)= X(J)*0.3048 
TP( J )= TP( J)/ 1. e 
TWP(J )= TWP( J)/1.8 

201 CONTINUE 

GAS= GAS ♦ 0.4535S24 
Q = 0 ♦ 1054.3503 
V( 1) = V( 1) ♦ 0.3048 
GASCHK= GASCHK * 0.4535924 
GSRATE= GSPATE ♦ 0.4535524 
WRITE (6,300) TIME, GAS, 0,V(l) 

WRITE (6,305) GA SCHK, G SRA TE , V ( NP ) 

WRITE (6,310) 

WRITE (6,311) (X( J),TWP(J),TP(J),J=l,N) 

210 CONTINUE 

00 211 J=1,N 
X ( J ) = X( J )/.3048 
TP( J) = TP( J ) * 1.8 
TWPU ) = TWP( J) * 1.8 

211 CONTINUE 

GAS = CAS/. 4535524 
0 = Q/1054.3503 
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V( u = V( D/.2C48 
GASCHK = GASCHK/. 4535924 
GSRATE = GSP ATE/. 4535924 
RETURN 
C 
C 

c 

1 FORf^AT (1HL/1HL»9H TIME * F6.2,9H SEC ONOS ,5 X , 17H GAS SUPPLIED = 

I F7.4,5H LRSf5X,17H HEAT TO WALL = F8.1.5H 8TU,5X, 14H INLET VE 

2L » F6.4, lOH FEET/SEC ) 

2 FORMAT (1HK,16X,17H CHECK CN GAS = F7.3t5H L8S,5X,13H GAS FLOW 

1= F7.4,9H LBS/SEC» 5X,15H OUTLET VEL = F6.4,10H FEcT/SEC) 

3 FORMAT (IFK.lOH X, FEET,2X»11H WALL T(R),11H GAS T(R),7X,10H 

1 X, FEET,2X, IIH WALL T(R),11H GAS T(RJ»7X,10H X, FEET,2X,11H 

2 WALL T(R),11H GAS T{R»/1X> 

4 format (F10.2.F 12.1»F 12.1 ,F 17.2,F12.l,F12.l.F 17.2»F12.1,F12.l) 

300 FORMAT I1HL/1FL,9H TIME = F6.1,9H SECONDS »5 X , 16H GAS SUPPLIED = 

IF7.4, SF KILOGRAM, IX, 17H HEAT TC WALL = E11.4,6H JOULE, IX, 14H INL 
2ET VEL = F6.4,10H METER/SECI 

305 FORMAT (1HK,16X,17H CHECK ON GAS = F7.3,4H KG,6X,13H GAS FLOW = 

1 F7.4,9H KGS/SEC , 5X, 15H OUTLET VEL = F6.4,11H METER/SECI 

310 FORMAT (1FK,10F X, METER, IX, 13H WALL T£MP-K,14H GAS TEMP-K,3X 
1, lOH X, METER, IX, 13H WALL TEMP-K,14H GAS T EMP- K, 3X , 10 H X, ME 
2TER,1X,12E WALL TEMP-K,14H GAS TEMP-K/1X> 

311 FORMAT <F10.2,F12.1,F 12.1,F17.2,F12.1,F12.1,F17.2,F12.1,F12.U 
END 
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Figure 3. - Logic diagram of expulsion program. 










Figure 4. - Logic diagram of ramp program. 
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